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APPLICATIONS OF THE CONJUGATE GRADIENT FFT METHOD IN 
SCATTERING AND RADIATION INCLUDING SIMULATIONS 
WITH IMPEDANCE BOUNDARY CONDITIONS 

Abstract 


The theoretical and computational aspects related to the application of the 
Conjugate Gradient FFT (CGFFT) method in computational electromagnetics 
are examined. The advantages of applying the CGFFT method to a class of large 
scale scattering and radiation problems are outlined. The main advantages of 
the method stem from its iterative nature which eliminates a need to form the 
system matrix (thus reducing the computer memory allocation requirements) and 
guarantees convergence to the true solution in a finite number of steps. Moreover, 
since the CGFFT algorithm is highly vectorizable, it can be efficiently implemented 
on supercomputers and multiprocessor machines. 

Results are presented for various radiators and scatterers including thin cylin- 
drical dipole antennas, thin conductive and resistive strips and plates, as well as 
dielectric cylinders. 

Solutions of integral equations derived on the basis of generalized impedance 
boundary conditions (GIBC) are also examined. These boundary conditions can 
be used to replace the profile of a material coating by an impedance sheet or insert, 
thus, eliminating the need to introduce unknown polarization currents within the 
volume of the layer. Moreover, by applying these surface boundary conditions, 
the difficulties associated with the calculation of the Green s function aie avoided. 
Impedance boundary conditions of up to the third order are employed and shown to 
be compatible with the basic CGFFT formulation, allowing an efficient simulation 
of large coated structures and filled cavity- backed apertures by further reducing 
the memory demand. For the purpose of validation of these simulations, a general 
full-wave analysis of two- and three-dimensional rectangular grooves and cavities 
is presented which will also serve as reference for future work. 
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CHAPTER I 


INTRODUCTION 


Despite its long life in classical electrodynamics, the study of “Radiation” and 
“Scattering” has enjoyed renewed interest in recent years, particularly in connection 
with improved antenna designs required by the technological progress in radio com- 
munication, advances in radar signature analysis and control, and more recently, the 
growing computing power offered by high speed computers. These have contributed 
to the emergence of Computational Electromagnetics (CEM), the numerical study of 
electromagnetic wave phenomena. 

When the operating frequency is such that the object’s linear dimensions are 
comparable to the wavelength, the available high frequency methods are no longer 
applicable and more accurate formulations must be adopted. Moreover, these asymp- 
totic techniques are predominantly suitable for conducting bodies of canonical ge- 
ometries and shapes. Therefore, they cannot be employed for simulating material 
bodies which constitute modern composite vehicles and structures. Furthermore, 
because of an increase in the complexity of the formulations and corresponding lim- 
itations on justifiable approximations, the need for consistent and stable numerical 
schemes arises to ensure convergence of the solutions under consideration. These 
restrictions, coupled with the limitations on available computer resources (memory 
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and speed), represent a challenge in the modeling of large scale problems. It is this 
class of problems that this study attempts to address. 

1.1 Motivation 

In radiation and scattering, we are interested in defining the electromagnetic fields 
in the presence of a source distribution. The key to the solution of any such problem 
is a knowledge of the induced current density on the surface or in the volume of 
the antenna or scatterer. Once this is found, the radiated or scattered fields can be 
computed via the standard radiation integrals. 

The induced volumetric current density J on the body of the scatterer or radiator 
satisfies an integral equation which may be expressed in functioned form as 

£' = m (l.i) 

where £' is a vector representing the impressed field and *4 is an integrodifferential 
operator (functional) relating the impressed fields to the induced current. Tradition- 
ally, equation (1.1) is solved directly by discretizing the unknown current density 
and forming a linear system of equations. Typically, such a discretization results 
in a square matrix demanding a memory storage of order C?(./V 2 ), where N is the 
number of unknown coefficients in the current density expansion. The system of 
equations is usually solved by standard matrix inversion methods such as Gaussian 
elimination or LU decomposition. However, the limitations on available comput- 
ing resources (including memory and processing time) associated with the numerical 
formulation of large systems, limit the range of applicability of such direct methods 
to relatively low frequencies. In addition, for large scale simulations, the memory 
demand of these methods results in prohibitive storage requirements and, thus, tra- 
ditional matrix inversion approaches are not attractive and alternative methods are 



3 


needed. 

To address this ne* d, iterative approaches have been used by researchers. In 
iterative methods, an initial solution for the current distribution is assumed, and 
this is improved through successive iterations. The process continues until a pre- 
assigned accuracy (tolerance) is reached. The main advantage of iterative methods 
is that the calculations can proceed without a need to generate the system matrix, 
because iterative methods often require only the multiplication of matrices with 
vectors. This reduces the memory requirement to a lower order O(N) and therefore 
renders iterative schemes suitable for large scale simulations. Furthermore, while 
a matrix inversion approach may fail to yield an accurate solution due to a large 
condition number of the matrix operators, an iterative method in such cases merely 
requires more iterations before reaching convergence. 

In this study, we will explore the application of an iterative scheme, namely the 
Conjugate Gradient Method (CGM), in the solution of systems of equations arising 
in scattering and radiation problems. From its introduction nearly forty years ago 
[1, 2], the CGM has been of considerable interest to mathematicians and engineers, 
primarily because, in theory, it ensures convergence for arbitrary initial estimates-a 
feature not shared by many of the iterative algorithms used in the past. 

The guaranteed convergence of the conjugate gradient method, as well as its effi- 
cient storage requirements as an iterative scheme, are prerequisites for its application 
to general configurations of interest. Another advantage of the CGM, however, stems 
from an interesting property shared by the integrodifferential operators encountered 
in most radiation and scattering problems. For these problems, A is a convolution 
operator involving the induced current density and the pertinent Green’s function. 
Thus, by employing the convolution theorem, the evaluation of the functional reduces 
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to simple algebraic operations on the Fourier transforms of the convolved quantities. 
This simplification often results in a notable improvement in the speed of the calcula- 
tions (Nlog N in contrast to N 2 ). In practice, the Fourier transforms are calculated 
efficiently via the fast Fourier transform (FFT) [3]. A CGM algorithm which in- 
corporates the FFT to carry out the convolution operations is often referred to as 
CGFFT method of solution [4], 

Other iterative methods utilizing the FFT algorithm have also been applied to 
a number of scattering problems [5]-[7]. However, these solution techniques usually 
suffer from two major deficiencies common to most iterative approaches: 1) conver- 
gence is not strictly guaranteed, and 2) convergence is often slow. The conjugate 
gradient method virtually eliminates the first problem because it guarantees mono- 
tonic convergence throughout the process. As for the second problem, the number of 
iterations required for the conjugate gradient method to yield a reasonable accuracy 
is often a fraction of the total number of unknowns. This depends primarily on the 
distribution of the dominant eigenvalues of the operator projected onto the system 
matrix. It has been argued convincingly [8] that the standard conjugate gradient 
method requires roughly twice as much computation time per solution as the Gaus- 
sian elimination, which is an 0(N 3 ) operation. However, the CGFFT is considerably 
faster since it requires only 4A^(1 -f log 2 N) operations per iteration-an overall 0(N 2 ) 
operation. 

The speed of the CGFFT method can be improved further by incorporating the 
subsectional expansion (basis) functions into the algorithm. In direct methods, it is 
well known that the use of appropriate expansion functions to represent the unknown 
current distribution plays an important role in the accuracy and convergence of the 
solutions. In fact, a large body of literature exists on various types of basis functions 
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and their implementation and performance in connection with direct approaches [9]- 
[11]. This is not the case with regard to the CGFFT method and, therefore, it is of 
interest to study the incorporation of these functions in the context of the CGFFT. 
It was found in the course of this research that such a treatment results in a drastic 
improvement (up to 100 percent) in the convergence rate of the method depending 
on the type of basis that is employed for the current expansion. 

Another area of interest addressed in this study is the incorporation of the gen- 
eralized impedance boundary conditions [12, 13] in the CGFFT method. General- 
ized Impedance Boundary Conditions (GIBC) are higher order boundary conditions 
which involve derivatives of the fields beyond the first. They have been found to 
be more effective than the traditional first order (standard) conditions in modeling 
thick dielectric coatings and layers [14]. The GIBCs have been successfully utilized 
in a number of analytical and asymptotic applications such as the Weiner- Hopf tech- 
nique and function theoretic approaches [14]. However, their utility in numerical 
methods has not been studied beyond the first order [15]. Applying these conditions 
on the surface of a dielectrically coated scatterer circumvents the need for sampling 
within the target’s volume and hence considerably reduces the number of unknowns 
required in the discretization of the problem. However, solution of these problems by 
direct methods is challenging due to the difficulty in handling higher order deriva- 
tives involved in the formulation. On the other hand, when the CGFFT method is 
employed, the derivatives may be carried out in the transform domain without much 
difficulty. These features make the formulation of such structures by generalized 
boundary conditions highly desirable and, therefore, a part of this study is devoted 
to the implementation and numerical study of GIBCs in connection with the CGFFT 


method. 
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Before closing this section, we remark that an inherent limitation of the itera- 
tive solution methods is their requirement that the solution process be repeated for 
each excitation. In direct approaches, on the other hand, once the system matrix is 
inverted, the solution for any excitation is virtually at hand. For this reason, itera- 
tive solution methods may be computationally intensive in those scattering problems 
where the target’s responses to a number of different excitations are of interest. This 
may be the case, for a example, when generating the backscatter pattern for an 
object, where the iterative solution must be repeated for each excitation. A partial 
remedy in this case is to use the results of the previous excitation as the starting 
point (initial guess) for the solution of the next excitation. For a certain class of large 
problems, however, the memory consideration may outweigh the possible disadvan- 
tages in speed associated with the multiple excitations. Therefore, in choosing an 
iterative method for these problems, the intensity of the computations and accuracy 
requirements as well as the merits of low memory allocations offered by such meth- 
ods must be carefully examined. In this regard, some vector and parallel processing 
features offered by modern computing facilities are important in reducing the CPU 
time in reaching convergence. At any rate, for problems involving a single excitation 
such as antenna radiation problems, the CGFFT is generally much faster than the 
general purpose matrix inversion techniques. 

1.2 Scope 

This dissertation is divided into three parts. Part One (Chapters II through IV) 
presents the CGFFT method as applied to radiation and scattering. In particular, 
Chapter II discusses the basic formulation and the incorporation of the subsectional 
expansion functions into the CGFFT method. Chapters III and IV present applica- 
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tions of the method to one and two dimensional problems classified according to the 
dimensionality of the employed Fourier transforms, respectively. These include 

• Radiation from thin wire dipoles, 

• Scattering by flat and cylindrical strips, 

• Radiation of cylindrical reflector antennas, 

• Radiation of dipoles in the presence of flat plates, and 

• Scattering by dielectric cylinders. 

For each application, the pertinent integral equations are derived and placed in a form 
suitable for a solution via the CGFFT method. Chapter V presents the generalized 
impedance boundary conditions and their incorporation into the CGFFT formulation 
for the simulation of two- and three-dimensional impedance sheets. 

Part Two presents a general study of a class of cavity structures and their anal- 
ysis using CGFFT in conjunction with the GIBCs. In particular, Chapter VI is a 
study of two-dimensional grooves of infinite extent, while Chapter VII presents a 
corresponding study of three-dimensional cavities recessed in perfectly conducting 
ground planes. 

Part Three discusses a vector- concurrent implementation of the CGFFT algo- 
rithm on supercomputers and multiprocessor machines. Chapter VIII presents re- 
sults from a numerical implementation of an optimized CGFFT algorithm, which 
further illustrate the potentials of the CGFFT in solving large electromagnetic scat- 
tering and radiation problems. 

In most cases, the accuracy of the solutions is confirmed by a comparison of 
the obtained results with available measured data or data obtained from alternative 
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solution techniques. Some of the presented results are out of the reach of direct 
solution techniques and can, thus, serve for validating future methodologies for large 
scale electromagnetic simulations. 

1.3 Basic Concepts 

When an object is exposed to electromagnetic fields, the scattered field U' is 
defined as the difference between the total field U T in the presence of the object and 
the incident field U‘ that would exist if the object were absent. That is 

U' = U T - U' (1.2) 

The fundamental laws governing the behavior of electromagnetic fields in space 
and time are Maxwell’s equations commonly expressed in differential form as 


V x E = 

dB 

dt 

(1.3) 

V x H = 

dD _ 
dt +J 

(1.4) 

V D = 

p 

(1.5) 

<1 

W 

II 

0 

(1.6) 


where 


E = Electric field intensity , volts/m 
H = Magnetic field intensity, amperes/m 
D = Electric Sux density (displacement), coulombs/m 2 
B = Magnetic flux density (induction), webers/m 2 
J = Electric current density, amperes /m 2 
P = Electric charge density, coulombs/m 3 
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and we shall use MKSC units throughout this study as indicated. The two curl 
equations (1.3) and (1.4) are Faraday’s induction law and the generalized Ampere’s 
circuit law, respectively, while the two divergence relations (1.5) and (1.6) are Gauss’ 
law for the electric and magnetic fields, respectively. The media interacting with the 
electromagnetic fields are characterized by the so called constitutive relations and can 
be classified according to their molecular structures and properties of their associated 
bound charge particles. For a sufficiently simple medium these relations are 


D 

= eE 

(1.7) 

B 

= pH 

(1.8) 

J 

= <7 E 

(1.9) 


where e, p and a denote the permittivity, permeability, and conductivity of the 
medium, respectively. 

In the presence of stationary material interfaces (the surface of a scatterer, say), 
an electromagnetic field satisfies the implicit(natural) boundary conditions 


n x (Ei - E 2 ) = 0 

(1.10) 

n x (Hi - H 2 ) = K 

(1.11) 

n • (Di - D 2 ) = p a 

(1.12) 

n • (Bi — B 2 ) = 0 

(1.13) 


where K and p a denote surface current and charge densities at the interface separating 
the two regions to which subscripts 1 and 2 correspond and n is the unit normal to 
the interface( usually taken to be outward with respect to the scatterer). The above 
boundary conditions are easily established on the assumption that the tangential 
components of D and B remain finite at the interface surface. It should be noted, 
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however, that the last two conditions are not independent of the first two for time- 
varying fields and are, therefore, redundant [16]. Moreover, if one medium is perfectly 
conducting {a — ► oo), no electric field exists in that medium as asserted by (1.9). 
Therefore, it follows from (1.10) that the tangential component of E is zero at the 
surface of a perfect conductor. 

In addition to (1.10)— (1.13), boundary conditions must be imposed at infinity 
to obtain unique solutions to the radiation problems. Physically, these radiation 
conditions require that solutions which represent outgoing waves traveling in a lossy 
medium vanish at infinity. 

Throughout this work we will consider harmonic time varying fields and adopt 
the time convention e ]Wt . Thus, the Maxwell equations (1.3) and (1.4) become 

V x E = — ju/B (1.14) 

V x H = jwD + J (1.15) 

with the explicit time dependence suppressed. The constitutive parameters under 
the time harmonic assumption are, in general, complex quantities. In particular, for 
a conducting medium, Ampere’s law (1.15) reads 

V x H = J e + J c + juT> (1.16) 

where J e represents the externally impressed current source and J c is the conduction 
current generated in the medium due to ohmic loss. Employing the constitutive 
relations (1.7) and (1.9), equation (1.16) can be rewritten 

V x H = J e + ju{t - ><t/u;)E (1.17) 

where, the quantity 


Cc = e — ja/u 


(1.18) 
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may be identified as the complex permittivity of the medium. 

The solution of a scattering problem consists of finding the solution of Maxwell 
equations which satisfy the boundary conditions at the surface of the scatterer and 
which displays the proper behavior at infinity. Typically, this is carried out by deriv- 
ing a suitable integral equation in terms of the unknown current density excited on 
the scatterer. Two popular integral equations for the time-harmonic electromagnetic 
fields are the electric field integral equation (EFIE) and the magnetic field integral 
equation (MFIE). The EFIE enforces the boundary condition on the tangential elec- 
tric field and can be used for both closed and open surfaces. The MFIE, on the other 
hand, enforces the boundary condition on the tangential components of the magnetic 
field and remains valid only for closed surfaces [17]. Since we are interested in both 
types of scatterers, the EFIE is developed and applied in this study. 

An equation for E may be obtained from Maxwell’s equations by eliminating H in 
(1.14) and (1.15) and using (1.7)-(1.9). Thus, assuming a homogeneous surrounding 
medium, we have for the scattered electric field 

V x V x E 8 - ifc 2 E 8 = -junJ (1.19) 

where k = u)y/fie is the wave number in the medium. This is known as the vector 
wave equation and the solution may be expressed as [18, 19] 

E 8 = -jkZ JJJ f (r; r') • J(r') dv' (1.20) 

where Z = yjpjt is the intrinsic impedance of the medium, f denotes the electric 
dyadic Green’s function in unbounded space given by 

r(r;r')=(l+pVv)G(r;r') 


( 1 . 21 ) 
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and G is the scalar Green’s function 


G( r; r') = 


-jk\r - r / | 
47r|r — r'l 


In the above, r and r' denote the field and source points, respectively and an explicit 


expression of T in Cartesian coordinates is 


a + JL & i & 

[ + k 2 dx 2) k 2 dxdy k 2 dxdz 


r= 1JL- 

k 2 dydx 


(1 + W> G <^'> 


(1.23) 


I 1 ^ 1 a 8 _ 1 _ 8» I 

\ k 2 dzd x k 2 dzdy + k 2 dz 2 / 

The magnetic field is then obtained from Faraday’s law (1.14). 

Equation (1.20) in conjunction with the appropriate boundary conditions on the 
tangential component of the total electric field gives the EFIE integral equation to 
be solved for the unknown current J. The specific form of this integral equation 
depends on the particular problem under study. Once the current distribution is 
evaluated from the integral equation, the scattered field throughout space may be 
calculated from the scattering integral (1.20). 

In radar applications, the target is completely characterized for the radar system 
by a quantity known as the radar cross section (RCS) or echo area denoted by <r(not 
to be confused with the conductivity & ). It is a measure of the reflective strength of 
the target and is mathematically defined as 


1 & 


a = lim 4nR 2 ^- 

R-+ oo P, 


(1.24) 


where P a is the power flux density of the scattered wave in a specified direction at a 
distance R from the scatterer, and P t is the power flux density of the incident plane 
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wave. The radar cross section is in general a function of frequency, polarization, 
and angle of incidence. When the incident and pertinent scattering directions are 
coincident but opposite in sense, (1.24) provides the monostatic or backscattering 
cross section. 

For two-dimensional targets which are infinite in extent along a given direction, 
the scattering parameter is referred to as the radar cross section per unit length or 
echo width and is defined as 

<r w = lim 2‘Kp 2 -£ • (1.25) 

p-*°o F, 

When the scatterer is long but finite in one dimension, the physical optics ap- 
proximation may be used to relate the three-dimensional radar cross section of the 
target to the associated two-dimensional echo width calculated on the assumption of 
infinite length. Hence, for plane wave Uumination normal to the long dimension of 
the scatterer, we have 

o 3d = 2{^) 2 <r u (1.26) 

A 


where £ denotes the length of the target in the long dimension. 
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Part I 

THE CONJUGATE GRADIENT 

FFT METHOD 



CHAPTER II 


THE CGFFT FORMULATION 


2.1 Introduction 

The integrodifferential equations considered in radiation and scattering have the 
general form 

E i (r) = ,(r)J(r) + j]j y f(|r - I'D • J(r')dv' (2.1) 

where E‘ denotes the excitation field, J is the unknown current density vector, f is 
the associated dyadic Green’s function, r and r' specify the observation and integra- 
tion points and t ) is the scalar function which depends on the electrical properties of 
the scatterer or radiator. 

In general, the above integral equation may be solved using direct methods such as 
the Method of Moments [20]. However, as the size of the problem increases, iterative 
techniques become more attractive for the solution of such equations. This is mainly 
because iterative methods avoid the process of matrix inversion which is subject to 
numerical instability for ill-conditioned matrices. Also, these schemes often involve 
only the multiplication of matrices with vectors and thus do not require an explicit 
storage of the system matrix. 


15 
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2.2 Description of the Conjugate Gradient Method 

The conjugate gradient method is a nonlinear semi-direct purely-iterative scheme. 
That is, assuming no truncation and roundoff errors, the exact solution is obtained 
in a finite number of steps depending on the number of independent eigenvalues 
of the operator matrix. This is achieved by applying the method to the normal 
equations obtained by premultiplying the system matrix by its adjoint. Moreover, 
the solution is improved at a monotonic rate throughout the process and convergence 
is guaranteed for a given number of unknowns and as the order of approximation 
is increased [21]. Convergence is accomplished via a systematic orthogonalization 
of the solution vector with respect to the residual vector defined as the difference 
between the left and right hand sides of the system at the end of each iteration. 
That is, for a system representing N unknowns, the solution vector is constructed 
from a set of N linearly independent (mutually conjugate) vectors orthogonal to the 
residual vectors. Since these also form a linearly independent set, the exact solution 
is obtained at the N - th iteration, but in general the solution can be constructed, 
rather accurately, with only a few of the orthogonal vectors that span the solution 
space. As a result, the desired tolerance is achieved in less than N iterations. 

The method starts out with an initial guess Jo and a corresponding residual error 
Ro. In each iteration, the residual vector is minimized not only along each local search 
direction but also over the entire span of search directions. To this end, the solution 
is expanded in terms of search vectors which would be generated by the modified 
Gram-Schrmdt orthogonalization scheme when applied to the sequence of residual 

vectors as the basis functions 1 . The set of search vectors {P n }, so constructed are 

x The choice of the TV-dimensional coordinate unit vectors as the basis functions would yield 
Gaussian elimination. 
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mutually A-orthogonal or conjugate (as opposed to orthogonal) 

<P„.4[Pj]>=0 , ijkj (2.2) 

The significance of this set of directions is as follows: for a quadratic function, 
successive line minimizations along a conjugate set of directions will achieve the 
minimum without the need to repeat minimization in any direction. Consequently, 
the minimum is achieved at the end of a finite number of steps. For nonquadratic 
functions, this guarantees quadratic convergence as the process goes on. 

A version of the conjugate gradient algorithm to be used herein is [4] 

R« = -4[J 0 ] - E* 

Po = -6-M a [Ro] 

Main Iteration Loop 

<n = PiFl 5 

Jn+l = Jn + LiP n 

R-n+1 — Rn + ^n-A[Pfi] (2-3) 

h 1 

n M“[Rn]l| 2 

Pn+1 = Pn — &n-4°[Rn+l] 

< 6 

Repeat If Necessary 

The norm and the adjoint operator are defined in terms of the inner product as 



||Uf =< u,u > 


( 2 . 4 ) 
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and 

< *4[U],V >=< U,«4°[V] > (2.5) 

It may be shown that [21] 

•4*W = l'(r) J(r) + JJJ v f(|r - r'|) • J(r') dv' (2.6) 

where * denotes the complex conjugate. 

2.3 Conjugate Gradient FFT Formalism 

The scattering integral in (2.1) is of convolution type and can, therefore, be 
evaluated in the spectral domain by invoking the convolution theorem. To describe 
this process, we must introduce the forward and inverse Fourier transformations. For 
one-dimensional functions, the Fourier transform pair is 

<){K) = j g(x)e~ jkxX dx (2.7) 

J —oo 

9 (*) = dk x (2.8) 

where k x is the spectral variable and we use the following symbolism to indicate the 


relationship among the transform pair 

g(x) <£> g(f x ) (2.9) 

Based on the above definitions, the convolution theorem is stated as [22]: 

f g(x')h(x - x')dx' = g(x) * h(x) g(k x )-h(k x ) (2.10) 

J — oo 

Similarly, the two dimensional Fourier transform pair is defined as 

g(k x ,k y )=[ [ g(x,y)e~ j{kxX+k ^dxdy (2.11) 

V) = 7 A 2 jQ f ^ 9( k x, k y )e i{kxX+k * v Uk x dk y (2.12) 
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with the convolution theorem expressed as 


— OO J —OC 


T 

g(x',y')h(x - x',y - y)dx = g(x,y) * hlx.yj -*=>■ g(k xl k v ) • h(k x ,k v ) 


(2.13) 


As usual, the transform of differentiated functions is given by 

<=► jk x g(k x ) (2.14) 

Using the Fourier transform notation, equation (2.1) can now be formally written 
as 

E* = t/J+ • J} (2.15) 

where T~ x denotes the inverse Fourier transform operator. Clearly, (2.15) avoids the 
generation of the square matrix corresponding to the operator A implying a storage 
requirement of O(N) as compared to 0(N 2 ) required with direct implementations. 
The solution of (2.15) via the CGM will be referred to as the CGFFT solution 
method. 

The Fourier transforms implied in (2.15) are, of course, continuous whereas in 
practice they will be replaced with discrete Fourier transforms (DFTs). It is, there- 
fore, necessary that an accurate relationship of the transforms in the discrete and 
continuous domains be established. Otherwise, a solution in one domain may not be 
representative of that in the other. Alternatively, excessive sampling may be required 
to represent the continuous system. 

Consider the finite-duration function g whose M consecutive sampled values cov- 
ering the entire domain of its definition are given by 


g n = g(x n ) 


x n = nAi 


n = 0, • • • , iV — 1 


(2.16) 
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The one-dimensional forward and inverse discrete Fourier transforms (DFT) of this 
sampled train are defined as [22] 


and 


9, = E 9«e~ i2m " N 

n = 0 


9,^ m ’ IN 

iV p= 0 


(2.17) 


( 2 . 18 ) 


where the consecutive spectral samples g v are separated by the spatial frequency 
interval A f x = l/(NAx). Similarly, the two-dimensional DFT pair is defined as 


and 


M-l N-l 

J W =EE g mn e-^r>IM +nq ,N) 

m=0 n=0 




(2.19) 


( 2 . 20 ) 


p =0 9=0 

For consistency, the transform of the differential operator (2.14) may be replaced 
by first approximating the operator by its discrete counterpart. For example, using 
a 3-point central difference scheme, we have 


dx [ ‘' - 


^£1 = ffjfi + 4?) -<y(j; - ^) 

.AxJ3 Ax 


( 2 . 21 ) 


whose Fourier transform is given by 

Tx ~ 9(k ’ ] 


_ e 


Ax 


( 2 . 22 ) 


or more compactly as 


da T 

f x <U jD x (k t )g(kx) 


(2.23) 


where 


At 

Dr = ) 


(2.24) 
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and sinc(x) = 8U j is known as the sampling function. It is seen that the transforms 
of the continuous and discrete derivatives (equations (2.14) and (2.23)) become equal 
as the spatial sampling interval becomes vanishingly small since 


lim sinc(x) = 1 

r — *0 


More accurate expressions may be derived by using higher order difference formulae. 
For example, employing the 5-point central difference scheme, we have 


dg , v 

aI (Xi) “ 


fA 9 


LAx 


_ + A£) - gjxi ± jjf) + gjxi £ ) - Mfi Z (2.25) 

5 1 2 Ax 


and the corresponding transform pair is given by 

dg T r 2 . ,. Ax. 1 . Ax.,^,, . 

4 r J^i[gSinc(A: I -^-) — -smc(3k x ^ )]<7(&x)‘ 


(2.26) 


2.4 Incorporation of Subsectional Expansion Functions 


In the Method of Moments, it is well known that the use of appropriate expansion 
(basis) functions to represent the unknown current distribution, plays an important 
role in the accuracy and convergence of the solutions. 

The employment of the subsectional basis functions to iterative methods involv- 
ing the FFT was initially proposed in connection with the Spectral Iterative Tech- 
nique (SIT) [23] and was shown to produce improvements in the rate of convergence. 
However, no quantitative conclusions were drawn because of convergence difficul- 
ties associated with the SIT. Here, a systematic study of the incorporation of basis 
functions into the CGFFT method will be considered. 

An assumption in the derivation of the DFT pair (2.17)-(2.18) is the validity of 
the integral approximation 


g{K) = f g(x)e ,k * x dx « Yldne. ]kxXn Ax 

J — oo 


(2.27) 
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implying that the integrand is constant over each sampling interval in (2.27). In 
other words, 

g(x)e~ }ksX = g rn + jg in = const. x n < x < x n+1 (2.28) 

where g rn = Re {<?„} and g tn = 5m {jr n }. A consequence of (2.28) is that g(x) is 
not constant over the interval and is, in fact, a function of both spatial and spectral 
variables. Thus, from a solution of (2.28) 

g T (x) = Re (< 7 (x)} = g rn cos(k x x) — </,„ sin(fc x x) (2.29) 

gi(x) = $Sm{g(x)} = g rn sm(k x x) + g in cos(k x x) (2.30) 


It has been observed [24] that the above dependence of the implied discrete rep- 
resentation of a given continuous function can play a major role in degrading the 
convergence rate of the CGFFT solution. It is, therefore, essential that some cor- 
rective procedure be taken and an obvious approach is to employ a higher order 
integration formula to replace (2.27). This was discussed in [24] but as can be ex- 
pected, it results in a slower DFT algorithm. An alternative [23, 25] is to expand 
g(x) in a sequence of subsectional expansion functions {/„} as 

N-l N - 1 

ff(x) = 52gnfn(x)= £ 9nf(x-X n ) (2.31) 

n=0 nsO 


where the second equality implies the invariance of the basis functions with respect 
to translation in the x-coordinate. Introducing the Dirac delta function 


S(x) = { 


1, x = 0 
0, else 

equation (2.31) may be rewritten as a convolution in the form 


(2.32) 


N - 1 

g(x) = f(x) *J^9n S{x - x n ) 

n= 0 


(2.33) 
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The Fourier transform of g(x) is thus given by 


9 = 1-9 


(2.34) 


where J is the Fourier transforms of the chosen basis function and g is the discrete 
Fourier transform of g as given by (2.18). Equation (2.34) establishes the relationship 
between the continuous and discrete Fourier transforms when subsectional expansion 
functions are employed. 

Customary forms of the basis function /(x) include the piecewise constant (PWC) 
and the overlapping piecewise sinusoidal (PWS) expansion functions given by 


P(x) = 


Q(x) = 


i 1*1 < 


Ax 


0 else 

sinffc 0 (Ax — |x|)] . I < ^ 

sin(fc 0 Ax) 111 S 


(2.35) 


(2.36) 


0 


else 


respectively where k 0 denotes the free space wave number. For these choices, we have 
the Fourier transforms 


— Ax. 

P(k x ) = Axsinc(fc x — ) 


(2.37) 


Q(k x ) = 


k 0 [cos(k x Ax) — cos(fc 0 Ax)] 


sin(fc 0 Ax) ( k * - k*) 
and we observe that for a sufficiently small sampling interval 


(2.38) 


lim P(k x ) = Ax 
a*— o ' 


(2.39) 


lim W*) 

&x->Q 


k 0 sm(k 0 Ax) — A: r sin(fc*Ax) 
cas(k 0 Ax)(kl — kl) 

Ax 


(2.40) 
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where use has been made of l’Hopital’s rule and the approximations 


sinx~x cosx~l ,|x|<Cl 


Therefore, for vanishingly small sub-intervals, the relation 


g(fr) czAx g 


(2.41) 


holds when the above expansion functions are employed. Apart from the multiplying 
constant, (2.41) is the transform of g(x) when f(x) = 6(x)-delta basis. Since the 
same result can also be derived via direct application of the rectangular rule of inte- 
gration (2.27) in the computation of the Fourier integral, (2.41) has been exclusively 
associated with the conventional application of the FFT algorithm despite the fact 
that it holds true for subsectional expansion functions as well. As will be shown 
later, the convergence of the CGFFT method is improved considerably if the more 
accurate expression (2.34) is used in the formulation instead of (2.41). 

In the case of a two dimensional current representations, an appropriate expansion 
is 

N-1M-1 

g(x, y) = f{x , y)*J2Yl - x n )6(y - y m ) (2.42) 

n=0 m=0 

where f(x,y) denotes the surface basis function and (2.34) still holds with the trans- 
forms interpreted as two-dimensional ones. 

Often, it is necessary that the basis function be chosen to have a different func- 
tioned dependence in the x and y directions. For example, when representing the 
currents on a thin plate a more suitable basis function is of the form 


f{x,y) = P{x)Q(y) or f(x,y) = Q(x)P{y) (2.43) 


having the corresponding Fourier transforms 
f = P(k x )Q(k y ) and 


}=Q(k z )P(k,) 


(2.44) 
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and P, Q, P, and Q are given in (2.35)-(2.36) and (2.37)-(2.38). Again, as Ax and 
Ay both go to zero, equations (2.44) reduce to 

g(k x , k y ) ~ ASg (2.45) 

where AS = Ax Ay is the incremental surface element. 

The above analysis enables one to incorporate the subsectional expansion func- 
tions into the CGFFT formulation. Using (2.34), equation (2.15) can now be written 
as 

E' = 77 • J + -J /} (2.46) 

Obviously, the transform f of the basis function needs to be computed only once and 
thus the computations per iteration implied by ( 2.15) and (2.46) are essentially the 
same. 

2.5 Numerical Considerations 

Equations (2.15) and (2.46) are valid only on the body of the scatterer. That is 
the domain of the Fourier transform is not infinite and for this reason, they cannot 
be solved directly for J in the spectral domain. To solve for J, (2.15) and (2.46) 
must be enforced on the scatterer in the spatial domain along with the sampling 
requirements and linearity of the corresponding discrete convolution [22]. In a dis- 
crete implementation of (2.46), the sampling intervals should be chosen so that the 
Nyquist criterion is satisfied in the spatial domain. Also the length of the FFT (re- 
ferred to as FFT pad) must be large enough to accommodate the spectral contents 
of the convolved quantities. That is, the truncation of the spectrum should cause 
minimal errors in the iteration process. 
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In general, the period M of the array to be transformed is chosen according to 
the relation 

N' = 2": N’>N Ny quut N' > 2N - 1 (2.47) 

where N is the number of unknown coefficients in the discretization of the current 
density and v is an integer. In practice, v is chosen to be the smallest integer 
satisfying the relation 

v ^ log 2 (2iV — 1) + e (2-48) 

where g is also an integer (usually unity) setting the order of the FFT pad. The array 
elements beyond the physical extent of the scatterer are set to zero before (forward) 
and after (inverse) transformation. 

In scattering computations, a usual practice in the implementation of (2.46) is 
to employ a sampling interval of at least 1/10 of a wavelength and an FFT length 
at least twice (order g = 1) that of the linear dimension of the scatterer in or- 
der to accommodate the spectral spreading due to the convolution. The sampling 
requirement is more serious for antenna problems where one is interested in an accu- 
rate evaluation of the surface fields for input impedance calculations. The FFT size 
should be chosen to minimize aliasing errors caused by the truncation of the Fourier 
transform of the Green’s function. However, as seen from ( 2.15) and ( 2.46), when 
minimizing aliasing, the entire quantity in the curly brackets must be considered. 
This involves the product of the transforms of the current with the Green’s function. 
When the current density is not expected to be associated with spatial singularities, 
its transform will be essentially bandlimited and an FFT length of order g = 1 should 
be adequate to represent the spectral content of the convolution without noticeable 
aliasing error. However, when the current density is associated with spatial singu- 
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laxities as in the case of E-polarized excitation for a thin conducting strip, aliasing 
is expected to cause substantial error unless corrective means are introduced. In 
general, to eliminate aliasing errors when employing the discrete Fourier transform, 
we must form periodic functions in the spatial and spectral domains [26] and this is 
the basis of the corrective procedure discussed later in the thesis. 

2.6 Summary 

A general overview of the conjugate gradient algorithm for solving electromag- 
netic scattering and radiation problems was presented. By introducing the Fourier 
transform pair and employing the convolution theorem, the electric field integral 
equation was placed in a form suitable for a solution via the conjugate gradient 
method. 

The incorporation of subsectional expansion functions into the CGFFT method 
was also discussed. A simple relationship between the continuous and discrete Fourier 
transforms of the unknown function was established in terms of the transform of 
the employed expansion function. The relationship holds for both one- and two- 
dimensional cases and may be considered as a generalization of a commonly used 
expression in the conventional application of FFT. The practical advantages of using 
the subsectional basis functions will be examined in the next two chapters. 

Finally, since the Fourier transforms involved in the calculations are computed by 
the fast Fourier transform, some numerical aspects of the method were also addressed. 



CHAPTER III 


RADIATION AND SCATTERING FROM 
WIRES AND STRIPS 


3.1 Introduction 

In this chapter the CGFFT method will be applied to the analysis of wire dipoles 
as well as flat and circular cylindrical strips. 

The radiation by a center-fed cylindrical wire dipole has been extensively studied 
with analytical approaches [27, 28] as well as traditional numerical techniques such as 
the method of moments [29, 30]. It is, thus, instructive to consider am application of 
the CGFFT solution method to this problem first. Two classic integral equations for 
the total current distribution over conducting wires are referred to as Pocklington’s 
integrodifferential equation and Hallen’s integral equation. The latter is usually 
restricted to the use of a delta-gap voltage source model at the feed of a wire antenna 
while the former is more general and is adaptable to other excitations. 

The scattering behavior of thin strips has also been studied in some detail in the 
last three decades. These include the scattering from conductive [31] and resistive 
strips [32, 33] as well as the analysis and synthesis of tapered strips [34, 35]. These 
studies have focused on flat strips. On the other hand, a numerical solution method 
for thin dielectric slabs of uniform thickness and arbitrary shape was given as early 


28 
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as 1965 [36, 37] by discretizing the slab and forming a linear system of equations to 
be solved via the Method of Moments. 

In the present study, the CGFFT will be applied for computing the scattering by 
the flat and circular strip problem. It is shown that for circular strips the convolu- 
tional form of the integral is preserved in terms of the angular parameter <j>. 

3.2 Radiation of a Thin Wire Dipole 


Consider a z-directed cylindrical dipole of length t and radius a < l radiating 
in free space. The electric field due to the excited current distribution K z over the 
antenna is given by the scattering integral (1.20). If the wire is thin, the current at 
the end faces is negligible and the radiated field in the cylindrical coordinates is then 
given by 




e -jk 0 R 

)- —adtfdz' 

1 4t rR * 


(3.1) 


where R is the distance between the observation point (p, <f>,z) and the source point 
(a, <£',*') 


R = \J(P -I- a 2 — 2 pa, cos(<£ — ^) + (z — z ') 2 


(3.2) 


and Z 0 and k a are the intrinsic impedance and wave number of the free space, re- 
spectively. For a (f> — symmetric method of feeding [38], the surface current density 
K z is azimuthally uniform and the total current is given by I z = 2iraK z . Moreover, 
since the radiated field is independent of (f>, we may set <f> = 0 for convenience. By 
enforcing the boundary condition 


Ei(a y z) + El(a,z) = 0, 


(3.3) 
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stating that the tangential field vanishes on the wire surface, we obtain the Pock- 
lington’s integral equation [39] 

E\{z) = jk 0 Z 0 (l + j-^j jT l z ( z ')G w (z - z')dz' (3.4) 

In the above, G w (z,z') is the Green’s function (also referred to as the exact kernel) 
given by 



(3-5) 


where 

R = ]j(z- 2 ') 2 + 4a 2 sin 2 | (3.6) 

The above integral equation may be simplified further for an electrically thin 
dipole(/; 0 a < 1). In this case, a total filamentary line-source may be assumed to 
flow along the center of the antenna along the z-axis. The angular integration over 
4> is avoided and G w is replaced by the reduced kernel 

e -jk 0 T 

G w (z - z') = — (3.7) 

4 irr ' 

where r is the distance from a field point (a,<f>,z) on the cylindrical surface to a 

source point (0, 0, z') on the z-axis 


r = \J(z - z') 2 + a 2 


(3.8) 


Comparing (3.4) with (1.1) we may identify the right hand side of (3.4) as A[I] 
whose adjoint is given by 



A form of (3.4) compatible with ( 2.15) is 

Fjz) = ^ JF-- {(*? - *J)G„(t.)/(t,)/} 


( 3 . 10 ) 
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where 


G w {k z ) 


~ k 2 0 )K 0 {a^k 2 z - k 2 0 ) exact kernel 
— K 0 (ayfk\ \ - k 2 0 ) reduced kernel 


(3-11) 


is the Fourier transform of the Green’s function in which I a and K 0 are the zeroth 
order modified Bessel functions of the first and second kind, respectively. Upon the 
specification of the excitation field E\ expression (3.10) is now suitable for a solution 
via the CGFFT method. 


3.2.1 Dipole Excitation Models 

Two excitation models commonly used in the analysis of the wire antennas, 
namely the voltage gap model and the magnetic frill model are considered here. 

Voltage Gap Model 

In this model we assume that the antenna is excited by a finite constant voltage K 
across its feed terminal gap giving rise to an impressed electric field which is entirely 
confined to the gap. Thus, the impressed field is expressed by 

E*' = zVJA \z\ < ^ (3.12) 

where A is the gap width. 

Magnetic Frill Model 

The delta gap model (3.12) does not account for the fringing fields present outside 
the gap region and, therefore, may not be accurate for near field and impedance 
calculations. Clearly, this situation becomes worse as the gap becomes wider. To 



Figure 3.1: The magnetic frill model for antenna excitation. 


include the effects of finite gap widths, the magnetic frill model was introduced [40]. 
This model is of practical importance specially in modeling of coaxial lines feeding 
monopoles on a ground plane (Figure 3.1). The feed terminal is replaced by an 
equivalent azimuthally directed magnetic current density that exists over an annular 
ring. The inner radius of the ring is chosen to be the same as the radius of the wire a, 
while the outer radius b is that of the coaxial cable feeding the monopole and whose 
characteristic impedance is 

ln(6/a) 

Zc = Zo-^ 1 (3.13) 


Assuming that the coaxial structure supports a purely TEM mode, its aperture field 
may be approximated by 


E* 


~ Vi 
^2pln(6/a) 


(3.14) 


and upon closing the aperture with a perfect conductor and invoking the equivalence 
principle in conjunction with image theory, we find that this field excitation can be 
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replaced by the equivalent magnetic current 


M = 2E' x n 


= -<f>- 


(3.15) 


P ln(6/a) 

The electric field generated by this source on the axis of the antenna is readily found 
to be [40] 

l e ~jkRi e -jkR 2 


E' z (0,z) = 


(: 

21n(6/a)' Ri 


R* 


■) 


where 


Ri = \/z 2 + a 2 R 2 = vz* + 5 2 


(3.16) 


3.2.2 Input Impedance 

Once the current distribution on the cylindrical body is known, the input impedance 
can be computed from 

1 r l 


z in = 


|/(0)l : 


J_E;(a,z’)r(z')dz' 


where E‘ x is the tangential surface field on the antenna given by 

E\{a,z) = -El 

Thus, the input impedance is given by 

1 r‘ 


z in = 


i/(o)i 


- j_P,{a,z')r{z')dz' 


(3.17) 


(3.18) 


(3.19) 


and for a voltage gap model, the above equation reduces to the well known Ohm’s 
law 


Z in = 


|/(0)l 


(3.20) 
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Figures 3.2 - 3.4 show results based on the above formulation along with compar- 
isons with data obtained by the method of moments (MoM). In particular, Figure 3.2 
and 3.3 show the convergence of the solutions as a function of sampling density using 
a magnetic frill model for the excitation fields and it is seen that the CGFFT and 
MoM solutions exhibit the same convergence characteristics. Also Figure 3.4 depicts 
the convergence of the input impedance (3.19) as a function of sample density and 
it is again observed that the CGFFT and MoM [41] solutions converge to the same 
result. 

The effect of incorporating various expansion functions is considered next. The 
current distribution on a 9A dipole based on a voltage gap excitation model is given 
in Figure 3.5 as predicted via a CGFFT or an MoM solution. Although all expansion 
functions considered give similar results, the employment of the piecewise sinusoidal 
basis functions (PWS) drastically improves the convergence of the CGFFT as seen 
from Figure 3.6. Typically, an estimated 100% improvement in the convergence rate 
of the CGFFT method was observed when employing the PWS expansion functions. 

Finally, Figure 3.7 shows the improvement in CPU time that can be attained on 
employing a CGFFT solution method versus a standard MoM solution. Clearly, the 
CPU time required for a CGFFT solution is a linear function of the system unknowns, 
whereas in the case of a MoM solution the dependence is quadratic. Also, shown in 
Figure 3.7 is the improved convergence attributed to the use of higher order basis 


functions. 



Figure 3.2: Numerical convergence of the linear current distribution for a 1A dipole 
with increasing sampling density evaluated by the CGFFT. Top to bot- 
tom: No. of samples = 15, 31, 63, 127; FFT pad order q = 2, 2, 2, 1; 
Magnetic frill excitation model. 
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Figure 3.3: Numerical convergence of the linear current distribution for a 1A dipole 
with increasing sampling density evaluated by the MoM. Top to bottom: 
No. of samples = 15, 31, 63, 127; Magnetic frill excitation model. 
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Dipole Input Impedonce; Frill Model 



Figure 3.4: Real and imaginary parts of the input impedance for a 1A dipole («/ A 
0.005) as a function of sampling density. 





Normalized Wire Current. |I(z)/H 0 | 
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Radiation of a Thin Wire Dipole 



Axial Position, z/\ 

Figure 3.5: Current magnitude for a 9A dipole (a = 0.005A) computed by the MoM 
and the CGFFT using different basis functions and a voltage gap model 
for the source (13 unknown/ A). 



Normalized Residual Error, As 
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Radiation of a Thin Wire Dipole 



Figure 3.6: CGFFT convenience patterns for the 9A dipole (13 unknowns/A). 
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Performance Evaluation of the OGFFT and MOM 


500 



Number of Unknowns, N 


Figure 3.7: A comparison of the CPU times required by the MoM and the CGFFT for 
the solution of the resonant dipole problem (CGFFT tolerance: 0.003). 
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3.3 Scattering from Flat Resistive Strips 

A thin conducting sheet or non-magnetic dielectric layer can be represented by 
a resistive sheet. In the case of a source- free dielectric layer having thickness r, we 
have from ( 1.17) 

V x H = jut c E (3.21) 

It is customary to characterize the layer by an equivalent electric current density as 

V x H = ju{t c — e 0 )E + ju>e 0 E 

= J e , + jve 0 E (3.22) 

where the equivalent current 

Je? — iwe 0 (e r — 1)E (3.23) 

is now assumed to radiate in free space. In the above, tr is the relative complex 
permittivity of the layer e c /eo- When the layer is electrically thin (kr <C 1), the 
normal component of the electric field inside the layer is negligible. The dielectric 
layer can therefore be replaced by a resistive sheet of surface current density 

K = limr[J„] tan (3.24) 

where 

[ J e?]tan = “ (» • J ««) » ( 3 - 25 ) 

is the transverse volumetric current flowing across the layer, (n is the upward unit 
normal to the layer). In view of (3.23) we may write 


E - (n • E) n = Z.K 


(3.26) 
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where Z, is the resistivity (in fi per unit squared) of the sheet 

= ~r r—T\ ( 3 - 27 ) 

jk 0 T{tr - 1 ) 

Therefore, a resistive sheet is an electric current sheet whose strength is proportional 
to the local tangential electric field. For a thin conducting sheet of conductivity a , 
(3.27) reduces to 

1 




<7T 


(3.28) 


Mathematically, the resistive sheet satisfies the boundary conditions [42] 

-in x n x (E + + E“) = Z.K 


n x (E + — E~) = 0 


(3.29) 


where E ± denotes the total field above and below the sheet. Using (3.29), integral 
equations may be derived for computing the current induced on the strips for a given 
excitation and in the following we consider their derivation and solution for each of 
the principal polarizations separately. 


3.3.1 Integral Equations 
E-Polarization 

Consider the E-Polarized wave 

JjW _ £gjM* “•*<>+ VMn^o) 

H‘ = — (x sin <t> 0 — y sin ^ o )y^ c , M IC “*° +v "“* ,<> ) 


(3.30) 

(3.31) 


incident on the resistive strip of resistivity Z s and width w coincident with the x-axis 
as shown in Figure 3.8. This excitation generates on the strip a z-directed current 
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Y 


Figure 3.8: Geometry of a strip illuminated by a plane wave. 

K z , giving rise to the scattered field 

E; = -jKZ. f n K,(x')G.(k 0 \x - x'\)ix‘ (3.32) 

J-w/2 

where G , is the two dimensional Green’s function given by 

G,(x - x') = jrHl 2 \k 0 \x - x'|) (3.33) 

V 

and Hj?) is the zeroth order Hankel function of the second kind. Imposing the con- 
dition (3.26) on the total tangential electric field over the strip, an integral equation 
for K t is obtained as 

L rw/2 

Y o< Jko* c~+ 0 = Kz ( x) + ^ Kg(x')HW (k a \x - x'|) dx‘ (3.34) 

4 J-w/2 

where r/ a = Z,!Z 0 is the normalized surface impedance of the strip. 
H-Polarization 

Consider now the H-polarized plane wave 

jji _ JgiMsco^o+vain^o) 

E‘ = Z 0 (x sin 4> 0 — v cos <j> 0 ) e 



jk 0 (x cos 4>o+y «in 4>o) 


(3.35) 

(3.36) 
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incident on the resistive strip. This excitation generates an x-directed current density 
responsible for the scattered field given by 

E‘ = -iKZ. (l + jf^ K,(x')G.(k.\x - x'\)dx' (3.37) 

Again, by imposing the resistive boundary condition, the integral equation satisfied 
by the current density K x is obtained as 

sin = r,.K,(x) + 7 fl + jT^ K x (x')U™(k,\x - *'|)<fx' 

(3.38) 

The far zone scattered fields at the cylindrical point (p, <f>) can be computed from 
the scattering integral using the large argument approximation of the Hankel function 

H™(k oP ) ~ * oP -* 00 (3-39) 

Upon using the approximations 

|p — p'\ ~ p — x r cos<j> ss p (3.40) 


for the phase and amplitude considerations, respectively, we have 

t~’ kop k 0 Z e 


El = 


V~P 


(3 - 4i) 


and 


E} = 


e _jfcp ,2LkZ fl 


1*7 

■e J 4 


\H-sin <f> r /2 A:x(x')e jfcoI ' c “^' 
V J—w/2 


(3.42) 


for E and H polarizations. The two-dimensional scattering echo width is defined as 


<12 


IF, 

a = lim 2irp , 

P-*°° ^| E *| 2 


(3.43) 
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and therefore 

. 2 

Z 0 r K t {x')e’ ko * ,a *+dx' (3.44) 

J-w/2 

. 2 

sin ^ [ K x {x')e ik ° x '™+dx' (3.45) 

J—w/2 

Typically, a solution of (3.34) and (3.38) can be accomplished numerically. How- 
ever, approximate analytical solutions exist for the perfectly conducting case if the 
strip is electrically very narrow or very wide. These solutions are based on the 
quasi-static and physical optics approximations of the pertinent integral equations, 
respectively. They may be used to find closed form expressions for the echo width 
of the strip. 

For a perfectly conducting strip Z t = 0 and the integral equations for the surface 
current densities are given by 

e i^co.^ = koZo r' 2 Kz ( x ’) H W( ko \ x _ x '\) dx ' (3.46) 

4 J-w/2 

and 

sin = Ml + jT* K.(x')H<?Xk,\x - x'\)ix' (3.47) 

3.3.2 Very Narrow Strips 

A general analysis of narrow strips and slots can be carried out analytically by 
employing certain quasi-static approximations to the integral equations developed in 
the previous section [43]. Since a similar analysis will be carried out in the study of 
narrow filled grooves in Chapter 6, we will present it here for completeness. 

When k Q w <C 1 in the integral equations (3.46) and (3.47), we may introduce the 
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small argument expansion for the Hankel function [44], 

.2 


H^\z) ~ 1 - j- In (Ip) + 0(z 2 , z 2 In z) 


(3.48) 


where In 7 = 1.78108 ... is Euler’s constant. Retaining only terms to 0(k a w ) in the 
Hankel function as well as the incident fields, we have 

C ln 11 - xVi ' (¥) + 4 C KAz ' )dx ' (3 - 49> 


for E-polarization and 

d 2 r> 2 

-w/2 

for H-polarization. Further, by introducing the change of variables 


d 2 /W 2 ; 

• 7 ^ J K x (x') In \x — x'\dx' = 2Trjk a sin <f> 0 


(3.50) 


f-2=. p-Sl 

w w 


(3.51) 


equations (3.49) and (3.50) respectively become 

/; - m' - ^ - [lo (if) + f ] /; (3,2) 


-£2 P KM') In — £'|d£' = jick 0 w sin <f> 0 


(3.53) 


To solve (3.52) and (3.53) we recall the following identities from the finite Hilbert 
transform theory [43, 45]: 

r 1 In jx — z'l , , 


f 1 In \x — x \ , , 

L Ji- X * dz = ~ lrh>2 [-l.i] 


(3,4) 


and 


d 2 /-i , 

y_j vl - x' 2 In |x - x'\dx' = ir x€[-l,l] (3.55) 

and since the right hand sides of (3.52) of (3.53) are independent of £ we deduce that 


KM) = 


Xe 


Z 0 JT=? 


Z a 4 1 




(3.56) 



47 


and 



(3.57) 


where \ e and \h are constants to be determined. By substituting (3.56)-(3.57) into 
(3.52)- (3.53) we readily obtain 


and 



(3.58) 


Xh = jk 0 w sin <f> 0 


(3.59) 


As expected (3.56) and (3.57) display the familiar edge behaviors at the terminations 
of the strip [46]. 

The scattering echo widths are computed from (3.44) and (3.45). However, in 
this case, we may use the approximation (3.40) for both amplitude and phase due 
to the small width of the strip. Thus, 


= 


O’A = 



(3.60) 


(3.61) 


and upon substituting for the currents, we obtain the simple expressions 



and 


<?h 


K 


1TW . . 

-g ~Xh sm <(> 



(3.62) 


(3.63) 


valid in the backseat ter direction. 
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3.3.3 Very Wide Strips 

For electrically wide strips, the local electric current may be assumed to be that 
corresponding to an infinitely wide strip. This is known as the physical optics ap- 
proximation and is expressed as 

K = 2n x H (3.64) 

or more explicitly, 

K z (x) = 2Y a sin <f> 0 e? k ° xco, ' t ’ 0 (3.65) 

for E-polarization and 

K g (x) = 2e? k ° xc '*+° (3.66) 

for H-polarization. 

The physical optics approximations may also be derived directly from the gov- 
erning integral equating (3.46) and (3.47) when the strip is assumed to be infinite in 
extent. Hence, we have 

e jkoX coe 4>o = jkoZo iim f w/2 Kz ( x ' )Gt ( x . x ' )dx ' (3.67) 

k Q w^oo J- w /2 

for E-polarization and 

sm<j> 0 e 3koXcr ** e = ^ Jto. £'>•(-') [(*1 + ~)} G.(^W (3.68) 

for H-polarization. The integrals on the right hand sides of (3.67) and (3.68) are 
equivalent to 

f K z (x')Cx,(x; x')dx' (3.69) 

J — OO 




G,(x; x')dx* 


—oo 


(3.70) 
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and are identified as convolutions in the infinite domain. Thus, upon invoking the 
convolution theorem and using the transform pair 

e jk a x cob4>o 2irS(k x - k 0 cos<f> o ) (3.71) 


we have, by taking the Fourier transforms of both sides of (3.67) and (3.68) 
2ir6(k x - k 0 cos <f> 0 ) = jk 0 Z 0 K z (k x )G a (k x ) 


(3.72) 


2jrsin <f> 0 8{k x — ko<x>s<j> 0 ) = T~K x (k x )(kg — k\)G t {k x ) 

Kq 


(3.73) 


Formally, the above equations can be solved algebraically for the transforms of the 
currents K z and K x to yield 

2tt 6(k x - fepcos^o) 


K z {k x ) = 


jk 0 Z 0 G a (k x ) 


and 


K t (k x ) = 


2 irk 0 sin <f> 0 S{k x — fcpcos^o) 


m - ki)G,(k x ) 

Taking the inverse transforms of both sides now gives 


«■> - 


p jk 0 x co* 4>o 


jk 0 Z o G,(k o cos<j> 0 ) 


W. n SQcz-kcco 

’ i J-°° <H - kl)G.(k.) 


(*i - 

k„ sin <f> 0 ei koXCt **° 


(3.74) 


(3.75) 


(3.76) 


(3.77) 


j(kl - kl)G t (k 0 cos<t> 0 ) 
where use was made of the properties of the 6 function. The Fourier transform of 
the Green’s function G a is given by (Appendix B) 

1 


G a (k x ) = 


2 


(3.78) 
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and when this is substituted in (3.76) and (3.77), we recover (3.65) and (3.66). 

To find the physical optics echo widths, we substitute (3.76) and (3.77) into (3.44) 
and (3.45). In the backscatter direction <f> 0 = <f > , we readily find 

a = k a w 2 sin <f>sinc 2 (k 0 w cos <f>) (3.79) 

for both polarizations. 

3.3.4 CGFFT Solution 

We will now consider the solution of (3.34) and (3.38) for arbitrary size strip 
via the CGFFT method. To do this, we must rewrite these equations in a form 
compatible with ( 2.46). The Fourier transform of G, is given by (3.78) 

G.(k x ) = - 

2 j 

and therefore (3.34) and (3.38) may be rewritten as 

y oejkoX = V ' Kz ( x) +jko jr-i {G.(k x )K z (k x )f(k s )} (3.80) 

and 

sin = r} t K x (x) + ■£?-' {(* 2 - kl)G t {k x )K x f{k x )} (3.81) 

respectively. These may now be solved via the CGFFT algorithm. 

Echo width patterns based on a CG solution of (3.80) and (3.81) axe compared 
with MoM data in Figures 3.9 and 3.10, respectively. The strip is 4A wide and has 
a non-uniform resistivity as shown. In practice, tapered resistive cards are often 
employed for radar cross section reduction and Figure 3.11 demonstrates an example 
of such a reduction in connection with a strip having a resistivity that is tapered 
parabolically as given in Figures 3.9 and 3.10. The choice of basis functions is again 




Normalized Resistivity, i) 











BlsUtlc Echowidth/X . dB 



Figure 3.10: H-polarization scattering results for a 4A parabolically tapered strip. 
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a factor in the convergence of the CG solution and similarly with the wire example, 
the sinusoidal basis functions were found to provide a substantial improvement in 
the convergence rate (almost 100 percent). This is shown in Figure 3.12 for the 
H-polarization. 

It should be noted that the expected behavior of the current density plays a 
major role in the choice of the FFT pad used in the calculations. This is related to 
the spectral content of the current as well as the singularity of the pertinent Green’s 
function. The field distributions over open conducting bodies and their singular 
behavior have been studied by several authors in order to establish such behavior in 
explicit numerical terms (see for example [47]). 

Consider the H-polarization case (TE Z ) first. In this case, the current density is 
not singular and-like the current density on the wire dipole studied in the previous 
section-it vanishes at the edges, rendering its transform essentially band-limited. 
Therefore, an FFT length of order 1 ( 0 = 1 ) should be adequate to satisfy the 
spectral spreading due to convolution without noticeable aliasing error. 

On the other hand, for the E-polarization incidence (TM Z ) the current density 
is singular at the edges and aliasing is expected to occur in the transform domain. 
This may cause substantial error unless high sampling rates are employed in the 
spatial and spectral domains to avoid aliasing. For example, the CGFFT solution 
for the perfectly conducting strip presented in Figure 3.9 for the E-polarization case 
required an FFT pad of order g = 3. This is, of course, undesirable because it will 
increase the memory deman o. and execution time per iteration. 

As mentioned in the previous section, employing the analytical transform of the 
Green’s function is valid if the inte al equation is defined on the entire real axis 
(infinite domain). This was the case for the wide strip in the limit as the width 



Backscatter Echowidth/X . dB Backscatter Echowldth/X . dB 
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Figure 3.11: Comparison of the backscatter echo widths of a 4A perfectly conducting 
and parabolically tapered strips. 



Normalized Residual Error, Xs 


55 


Scattering from a Conducting Strip 



Figure 3.12: Convergence patterns for the 4A strip illuminated by an H-polarized 
plane wave using 20 unknowns/ A. 
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was taken to infinity. For finite strips, however, this approach gives an approximate 
solution which improves by extending the size of the FFT pad to include higher 
spectral components. For a given strip, the degree of improvement achieved depends 
on the sampling density and the polarization of incident field. 

To overcome this difficulty, an alternative is to discretize the integral in (3.34) 
before proceeding with its computation via the discrete Fourier transform. That is, 

assuming a pulse basis expansion for the current density, 

N-l A 

Kz{x) = K x {x n )P{x -x n ), x„ = n A + — (3.82) 

n=0 * 

we substitute (3.82) into the integral equation (3.34) and enforce it at discrete 
points(point-matching) x m = m A + m = 0 , ...,N - 1. This yields a linear 
system of equations for the solution of the current density. In particular, we have 

rw /2 ru //2 I **- 1 

/ K,( x ')HfHK\x m -x‘\)ix' = / //< 1 >(Jfc„|x„-x'|)<ix' 

J-w/2 J- w /2 [ n=Q J 

(3.83) 


and by interchanging the order of summation and integration, 

1 N — l 

E *.(*.) / ^ (fcl*. - A)ix’ = E - X.) (3.84) 

n =0 J x n ~~ 2 n= 0 

In the above, T(x m — x n ) = T mn are the mutual admittance elements given by 


j2 /.At A 




T 

A mn — 


(3.85) 


^Atfi ! ’(*„|x m -x„|), 


n / m 


in which In 7 is the Euler’s constant. Since T is not singular anywhere, the evaluation 
of the convolution integral may now be carried out without aliasing errors via the 
discrete Fourier transform as 


K t (x')H ^ ] 'Alxm - x'\)dx' = DFT - 1 [k z f } 


(3.86) 
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where DFT -1 denotes the inverse discrete Fourier transform and T is the discrete 

transform of the sample train T on, n = — ( N — 1), • • • , N — 1. 

Expression (3.86) renders the evaluation of the convolution relatively insensitive 
to the length of the FFT provided the convolution requirement is satisfied. As 
illustrated in Figure 3.13 for the case of normal incidence on a perfectly conducting 
strip one wavelength wide, the predicted current distribution agrees with the MoM 
result when (3.86) is employed in the CGFFT algorithm with an FFT size just twice 
the length of the strip (FFT pad of order g = !)• 111 contrast, when employing 
the sampled continuous analytical transform for the evaluation of the convolution 
integral, the resulting current distribution remains in disagreement with the MoM 
solution unless at least an FFT pad of order g = 3 (four times the size of the strip) 
is used. The corresponding comparison of the bistatic scattering patterns is shown 
in Figure 3.14 and the same observations again apply. 

Hereon, the solution of (3.86) via the conjugate gradient method will be referred 
to as the CGDFT method and the corresponding method of solution based on (3.80) 
or (3.81) will be referred to as the CGFT method. 
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Scattering from a Conducting Strip 

w=1.0X, E-Pol., Normal Incidence 



x/\ 


Figure 3.13: Comparison of the current distribution on a 1A wide perfectly conduct- 
ing strip illuminated by a plane wave (E-pol, <f> 0 = 0) as computed by 
various methods. 
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Scattering from a Conducting Strip 

w=1.0X, E-Pol., Normal Incidence 



Figure 3.14: Comparison of the bistatic echo width of a 1A wide perfectly conducting 
strip illuminated by a plane wave (E-pol, <f> 0 = 0) as computed by 
various methods. 
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Figure 3.15: Geometry of an infinitely long curved strip illuminated by a plane wave. 

3.4 Scattering from Cylindrical Strips 

Consider a thin cylindrical shell of resistivity 17 illuminated by a plane wave E* of 
wave number k 0 and polarization angle ip (Figure 3.15). The incident field is given 

by 

E*(p) = (^sin0 + zcostp)Z 0 e~jk<>(ki ■ p) (3.87) 

where ip = 0 corresponds to E-polarization (TM Z ) while ip = — tt /2 corresponds to 
H-polarization (TE Z ). 

The scattered electric field due to the excited surface current K on the shell is 
expressed by the line integral 

E» = - jk 0 Z 0 jf K(p') • f(p; p')dl' (3.88) 

where F denotes the electric dyadic Green’s function in unbounded space given by 

IW)= 



(3.89) 
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In the above, G t is the two-dimensional free space Green’s function given by 


4; 

The explicit form of f 1 in cylindrical coordinates is 


/ Id 2 11 

v 1 • L 1 a -2 / 1.2 


r = 


kid p*’ kl p dpd<t> 


1 1 d 2 


(1 + 


i d 2 


k 2 pd<j>dp v (k 0 p ) 2 d<t > 2 


) o 


0 


0 


(3.90) 


G s (p;p') (3.91) 


\ w / 

The total tangential electric field on the strip satisfies the resistive boundary 

condition (3.26) 


E t - (n • E t ) n = |E* + E*| = Z t K 


which upon substitution of (3.88) yields the desired integral equation 

- n x n x E i (p) = Z.(p)K(p) + jk 0 Z 0 j K (p f ) • T(p; p')dl' (3.92) 

to be solved for the unknown current distribution. In view of (3.90) this represents 
a convolutional integral equation in K. 

Let us now consider a solution of (3.92) for the special case of a circular strip of 
radius a. Referring to Figure 3.16 and defining the phase reference at the origin, we 
may write 


dl 

Ip - Pi 


ad<f> 


a [(cos 4> — cos + (sin <f> — sin <f>') 2 


2asin(| ^ |), 


1 

2 


(3.93) 
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Figure 3.16: Geometry of a circularly curved strip. 

and 

ki- p = p 0 -a cos (<f> - <(> o) (3.94) 

Also, since there is no variation in p and the strip is infinitesimally thin, (3.91) 
reduces to 

^ = 4 J { \ l + (koa) 2 tt + zz } H l 2) [ 2fc 0 asm(— - - -) ] (3.95) 

The integral in (3.92) can thus be expressed as a convolution in <f> yielding 

+ z - 4>o) = v .(<t>)K(<f>)+ jk 0 jK(t') ■?(<(>•, <t>')ad<t>' (3.96) 

where r) t is the normalized surface resistance of the strip and use has been made of 
(3.94). For TEz incidence this becomes 

E% M = + “J" 1 + ] J c [ 2k o“ s in( — 2 ~ - ) ] d<f>' 

(3.97) 

while for TMz incidence we have 



E'.W = 1 1.WK.W + [ 2t„a 8 in(Ji^il) ] df, 


(3.98) 
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Clearly, both (3.97) and (3.98) are amenable to a CGFFT solution. 

It is noted that if the radius of the strip is large compared to its width, we may 
modify the approximations (3.93) and (3.94) as 

lim | p — p'\ = a\<f> — 4>'\ — I* — x'\ 

a— +oo 

lim acos(<^> — <f>$) — a sin <j>o + a(7r/2 <^)cos<^o 

a— *-oo 

= as\n<f>o + xcos<f>o (3.99) 


and the formulation reduces to that of scattering from a flat strip (see (3.38) and 
(3.34)) 

(_* sin * sin t + 2 cos *)«>'*•<« sin * + 1 cos *>) = 

/ w/2 _ 

K(*0-r(*;*0 rf * 

- v //2 

where K(x) = xK x (x) + zK z (x) and f is now given by 


(3.100) 


f(x;*') = 4j 


1 d 2 

(i + — or;)* 2 + " 


H?\k 0 \x - x'|) 


(3.101) 


k*dx 2> 

To solve (3.96) the current density is expanded in terms of a subsectional surface 
basis function ^ as 

AT— 1 

(3.102) 


KM =i)K.. *.(*) 

n =0 


where 


$ n ($ = $(<£ - < t > n ) 


(3.103) 


and upon using the piecewise constant basis function, 


(3.104) 


with 


P{4) = 


l 

o 


\< f >/ 2 < ^ < A < f >/2 


else 


(3.105) 
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Substituting for the current expansion in the integral on the right hand side of (3.96) 
and interchanging the order of summation and integration, gives 


£k,. / *.(#•) -r(*#w 

Tl JC* 


(3.106) 


Introducing (3.106) into (3.96) and applying Galerkin’s method (Appendix A) yields 
the system of equations 


V m = A<t>r] tm K m + ^ £S ron • K n 


(3.107) 


where 


(3.108) 


V m = (4>sinxl> + z cos VO / ejkoCo S (<t> - 
The dyadic function 

H = (3.109) 

is a discrete kernel whose elements are given by 

r ^,+A/2 i d 2 ] L. . .\4>-<t>'\ 


tft = r f & a— (“) 

J<t>m-A/ 2 2 {k 0 a) 2 d<t> 2 J [ v 2 7 


f“n = r + 7T +A, >> 

- 4>Tk — i 


-A/2 


2Ar 0 a sin(— — — ) 


d4>'d4> 


d<j>'d<t> 

(3.110) 

(3.111) 


It is noted that both £ ^ and <f have integrable singularities corresponding to the 
self-cell interaction which can be approximated analytically. In particular, 


>nn — 


(k 0 a) 


r) 


v^FJb 0 aA^J 2 > - 2 j 


(3.112) 


and a similar expression may be obtained for The remaining terms (n ^ m) 
may be evaluated numerically. 

Applying the discrete convolution theorem in (3.107) now gives 


V m = A<t>r) m K m + ^ DFT -^I . K} 


(3.113) 
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which is in a form suitable for solution via the conjugate gradient method. 

Once the surface current density is evaluated, the scattered field can be computed 
using (3.88) as 

E*(*) = -jZ 0 k 0 a jK(<f>') • f(|* - 4>W 

Specializing this to the far-field, we find the scattering echo widths for the two 
principal polarizations to be 

= - I a / ~ 2 (3.114) 

4 I Jc 

for TM Z polarization and 

<7 h (4>) = ^ asin^o jKM’)J k ° acos (*- | 2 (3.115) 


for TE Z polarization. 

Sample calculations are now presented for circularly curved strips using the above 
formulation. Figure 3.18 shows the bistatic scattering patterns for a 2 A fiat strip as 
it is uniformly bent to form a closed circular cylinder keeping its width (perimeter) 
constant. The strip is positioned symmetrically around the y axis and illuminated by 
a TM Z plane wave incident at 90 degrees. It is noted that as the curvature k increases 
from zero (flat strip), the main (specular) lobe drops and eventually disappears in 
the limit when the complete cylinder is achieved. The numerical result for the closed 
cylinder is in agreement with the classical eigen-function solution [17] 


/ n 2 ^ 2 Jn{k 0 a) , 2 ,o n, 

; S i <3 - 111 

where J~ is the Bessel function of order n and 6„ is the Kronecker delta function 


| 1 P = 9 
1 0 9 


(3.117) 


This is illustrated in Figure 3.19 for the same strip. 
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Figure 3.17: A 2A wide conducting strip as it is uniformly bent to form a hollow 
cylindrical tube, k is the curvature of the strip and 6 is the polar angle 
subtended by the strip. 

3.5 Radiation by Cylindrical Reflector Antennas 

Consider the circular cylindrical reflector shown in Figure 3.20 ill umin ated by the 
line source 

E, = (3.118) 

The total electric field E T is evaluated in the far zon e(k a p » 1) as 

+ (3.119) 



67 


Bistatic Echowidth of a Circular Conducting Shell 
w=2A, Normal Incidence ^ 0 =90* 



Observation Angle <p % deg. 


Figure 3.18: The bistatic echo width of the strip in Figure 3.17 illuminated at normal 
incidence. 

with the normalized radiation pattern of the reflector antenna given by 

FM = 1 + / ■K-,(*V i - oa>s (* - *'W' f (3.120) 

i e Jc I 

Figure 3.21 shows the radiation pattern of an infinite electric line source in the 
presence of a 2ttA cylindrical resistive strip (a = 4A/3). The line source is positioned 
at the center of the strip and radiates through a right angle slit. As expected, the 
nonzero resistivity reduces the directivity of the reflector. 
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Bistatic Echowidth of a Circular Conducting Cylinder 
*r=2A, Normal Incidence <p o = 90 *, E-Pol. 



Figure 3.19: A comparison of the computed bistatic echo width of a circular con- 
ducting cylinder with the 20-term eigen-function solution. 
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Figure 3.20: Geometry of a cylindrical reflector antenna with a 90 degree circular 
slit excited by an infinite electric line source. 

3.6 Summary 

Scattering and radiation from thin wires and strips were formulated using a stan- 
dard integral equation approach. The convolutional integral equation was uniformly 
discretized allowing the implementation of the fast Fourier transform for carrying out 
the calculations. For the antenna problem, a larger sampling density was required 
to yield an accurate evaluation of the input impedance. 

Two formulations for a conjugate gradient solution of the scattering by resistive 
strips were presented. The first formulation, namely, the CGFT method employed 
the sampled continuous transform of the Green’s function for the evaluation of the 
convolution integrals. The other formulation, called CGDFT, employed finite dura- 
tion discrete Fourier transforms for the evaluation of the same integrals. This was 
found to provide a more accurate as well as a more efficient simulation since it elim- 
inated all aliasing errors. Notably, the system solved by the CGDFT method is the 
same as that generated by the standard moment method procedure. 


Normalized Pattern F B (<p), dB 








71 


It should be noted that (3.38) and (3.81) are also applicable for computing the 
scattering by an impedance insert of width w. This simply requires the replacement 
of Z a by the impedance of the insert and changing the polarization of incident field. 
The resulting echo width is then twice that given for the resistive strip to account 
for the presence of the ground plane. It should also be noted that the formulation 
discussed in connection with the thin strips is equally applicable to circular slabs 
of finite thickness by introducing equivalent volumetric currents instead of surface 


electric currents. 



CHAPTER IV 


RADIATION AND SCATTERING FROM 
PLATES AND CYLINDERS 


4.1 Introduction 

Planar and cylindrical structures constitute simple but nevertheless important 
components in man-made structures. Simulation of electromagnetic scattering from 
these targets is of academic interest as well as practical value in computational elec- 
tromagnetics. Understanding the electromagnetic scattering behavior of these struc- 
tures is also important in modeling more complex targets as well as in radar detection 
and cross section reduction. Although plates and cylinders have been the subject 
of intense study in a wide range of frequencies, their numerical analysis have been 
limited to the low frequency region, primarily due to computational limitations of 
the traditional direct methods. In particular, experience with various numerical and 
asymptotic methods of solution as well as comparison with measured data reveals 
that there is a serious difficulty in accurately predicting the scattering behavior of 
plates at grazing incidences where the edge currents and corner diffraction effects are 
significant. 

In this chapter, we first develop the necessary integral equations which are then 
transformed to a suitable form for a solution via the CGFFT method. Two ap- 
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proaches will be employed in the application of the method. The first implementa- 
tion, previously referred to as the CGFT method in connection with the strip anal- 
ysis, employs the sampled continuous Fourier transform of the free space Green’s 
function for the evaluation of the pertinent convolution integrals. This approach 
assumes an infinite spatial domain in the definition of the Green s function. Thus, as 
far as the Green’s function is concerned, the finiteness of the target s physical extent 
is not accounted for and unless a large FFT ‘pad’ with extended zero elements is 
used, the method suffers from aliasing errors. A pad at least 3 times the size of 
the target in each dimension is often needed to obtain acceptable results at oblique 
and close to grazing incidences [48]. To alleviate this difficulty, another approach, 
previously referred to as the CGDFT method in connection with the strip problem, 
will be employed where the pertinent integral equation is first cast into a discrete 
form before the application of the convolution theorem to evaluate the integrals. As 
observed in the case of the strip, this eliminates all aliasing errors, except perhaps 
those attributed to a possible under-sampling of the current density. 

Below, we discuss both of the above formulations for the solution of integral 
equations arising in the computation of the scattering by resistive plates and dielectric 
cylinders of arbitrary shapes and cross sections. The accuracy and efficiency of 
these formulations are then examined by a comparison with measured data and data 
generated by alternative techniques. 

4.2 Scattering from Thin Plates 

Consider a thin inhomogeneous plate of resistivity Z t illuminated by an incident 
field Ej and we are ini : ted in evaluating the scattered field from the plate. 

The scattered field due to the excited surface current density, K on the plate is 
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given by the surface integral 


E*(r) = —jk 0 Z 0 J' K(r') • f(r; r')*' (4.1) 

where T denotes the electric dyadic Green’s function in unbounded space given by 

f(nO=(l+iw)c F (r;0 (4.2) 


with 


G>(r;r') = 


o-jK |r - r'| 
4x|r — r'j 


The explicit form of f is now given by 

/ /t Id 2 , 


r = 


l d 2 


kldx 1 ’ k\dxdy 


1 d 2 


/, 1 d 2 , 

(! + TiaTa) 


G P ( r;r') 


(4.3) 


(4.4) 


V kldydx v kldy*' ) 

The total tangential electric field on the plate satisfies the resistive boundary 
condition (3.26) and the desired integral equation for the unknown current density 


is 

PwU - Z.MKM + A* l km ■ rte r')*' (4.5) 

in which r and r' denote the field and source points on the surface of the plate. 
Expanding the current density in terms of & subsectional surface basis function 
we write 

M-l N-\ 

K(x,y) = ^ £ K mn * (4-6) 

m=0 n= 0 

where 


*mn(*,y) = ^(x -x m ,y- y n ) 


(4.7) 


and 


®(x,y) = xxrf> x (x,y) + yyif> v (x,y) 


(4.8) 


in which an d t/> v are the expansion functions in the x and y directions, respectively. 
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4.2.1 Formulation Using Continuous Transforms 

Through application of the convolution theorem, the continuous transform of 
J(x,y) as given in (4.6) cam be written as 

J=jl (4.9) 


where J = xJ x (x,y) + yJ v {x,y ) denotes the two dimensional discrete Fourier trans- 
form of the train J m n defined in ( 2.19). Also, J(&x, ky) = xJ x {k x , k y ) -I- yJ v {k x , k y ) is 
the continuous Fourier transform of J defined in ( 2.11) and t denotes the continuous 
transform of the basis function. 

By invoking the relation ( 2.14), the continuous transform of the free space dyadic 
Green’s function can also be written as 


G p (k x ,k y ) (4.10) 

V K / 

where G p is the transform of the Green’s function given by (Appendix B) 


*1 V 


f = 


L L k 2 

rvj Ay . 

-an 1 • n 


G p (k x , fc v ) = 


(4.11) 


2 jjkl -kl-k* 

Equations (4.10) and (4.11) constitute analytical expressions for the Fourier trans- 
form of the free space dyadic Green’s function. Substituting these into (4.5) and 
testing the resulting equation at discrete points (point-matching), yields the system 


Etf = + jk 0 Z o r~ x {{T • #) • J} (4.12) 


where the subscript ij denotes the value of the quantity at the test point (xj,y : ) on 
the plate. It should be noted, though, that in performing the Fourier transformation 
implied by (4.12), an FFT pad at least twice the size of the plate in each dimension 
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must be employed. In general, however, a much larger pad is required when the 
analytical transform of the Green’s function is used. Also care must be applied when 
implementing (4.12) to avoid sampling at the singularity of the transform of the 
Green’s function as given in (4.11). 


4.2.2 Formulation Using Discrete Transforms 


In this formulation, the integral equation (4.5) is first discretized leading to ex- 
pressions that can be identified as finite domain discrete convolutions. These can 
then be evaluated via application of the discrete convolution theorem which is in- 
herently cyclic, thus, avoiding aliasing errors. To cast the integral equation (4.5) in 
discrete form we first employ (4.6) to rewrite the right hand side integral as 


L 


j-M-1 TV— 1 

E X) J mn ' $mn(*,y) 

m=0 n=0 


r (x,y;x',y')ds f 


(4.13) 


which, upon interchanging the order of summation and integration, may be written 
as 

M~1 T V-1 

(4.14) 


H ■f(x,y ) x',y')ds' 

71=0 n=0 •'*' 


Introducing (4.14) into (4.5) and satisfying the resulting equation at a discrete set 
of points (point matching) yields the system of equations 


®*»i - ZiijKy + j k D Z 0 ^ a,,- • K r 


(4.15) 


The dyadic function 

^ £ xx 

£V X £VV 

is a discrete kernel whose elements are given by 


(4.16) 


(ff = (* + p Gpfoyjx'yhM*' -x m ,y' -y n )ds' 
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«i = L m , Gp(l ’ v; x ’’ * ,) ' fc,( *’ - Xm ’ y '~ v " )ds ' 


(4.17) 


{”” = (l + p^5 )/ G r (x,y,x',y')^ v (x'-x m ,y'-y„)ds' 


where cr mn is the incremented surface element corresponding to the mnth cell on the 
plate and all expressions are evaluated at (x, y) = (x<, y,) upon differentiation. Obvi- 
ously, the convolutional nature of the operation is preserved once the above functions 
are evaluated at the appropriate field points. Applying the discrete convolution the- 
orem in (4.15) now gives 


E a = + jkZ 0 DFT "Ml • K} (4.18) 


where 5 denotes the discrete Fourier transform of 2. 

To calculate the elements of 2, the partial derivatives may be carried out by finite 
difference formulae. In particular, using a 3-point central difference scheme ( 2.23 ), 
we find that 


§ = 1 


1 kl-D 2 x —D x D y 


(4.19) 


V -D y D x kl-Dl) 
where l is the discrete Fourier transform of the sequence(assuming piecewise constant 
basis functions) 


£mn — / 


,-jkos/x 2 + y 2 
±Ty/x* + y 2 


■da 


(4.20) 


and 


D x 

Dy 


k x sinc(k x ^-) 

fc v sinc(fc v -^) 


(4.21) 

(4.22) 
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as given in ( 2.24). It is noted that £ has an integrable singularity when x m =y n =0 
corresponding to the self-cell interaction. This term can be evaluated analytically 
using one of the following approximations: 


Approximate integration : From [49] 


+ Ay In tan ( j + i tan"' - jk „ 


Ax Ay \ 


1 


Taylor series expansion : Expanding the integrand of (4.20) as [50] 


,-jk 0 R 

R 


“K 


1 — jk 0 R + 


UMl UMl 


) 


£oo can then be expressed as 


foo a -L(/, - jkj, - S/ 3 + 


2 ' ~ 


— * 




where 


h 

h 

h 

u 


J f = xln(y + R) + yln(x + R) 

J J ds = Ax Ay 

J f Rda = + y ln(y -f R) + y ln(x + R) 


Circular disk approximation: 


rr 0 e ~jk 0 r 

l L -rr^rdvi^, 


4?rr 


= 


(4.23) 


(4.24) 


(4.25) 


(4.26) 


where 


r 0 = \J Ax Ay/ ic 


(4.27) 



79 


Figure 4.1 shows a comparison of these expressions for £ 00 for square cells (Ai = 
Ay = A) of different sizes. As seen, they all give values that are essentially indis- 
tinguishable for A < 0.1 A. The remaining terms £ mn are evaluated numerically via 
Gaussian quadrature integration. Using the above formulations, computations were 
performed for a variety of plate sizes and shapes under two different excitations, 
namely, plane wave excitation and Hertzian dipole excitation. 

4.2.3 Plane Wave Scattering 

Plates have been of considerable interest in plane wave scattering because they 
often represent building blocks in the simulation of more complex configurations of 
practical interest. An understanding of their scattering characteristics can, there- 
fore, provide insightful information for design applications. In this case, simple high 
frequency formulae are usually more suitable, but unfortunately, available expres- 
sions have not been found to yield accurate results. On the other hand, numerical 
simulations demand an excessive storage requirement making the CGFFT solution 
method attractive for such simulations. 

Consider the plate in Figure 4.2 illuminated by a plane wave 

E 1 = E 0 e-i k °( ki ' r ) (4.28) 

H 1 = ki x E‘ (4.29) 

Zo 

where Z a and k 0 are the free space intrinsic impedance and wave number, respectively. 
In the above, fc, is the unit propagation vector 

% = — [sin 0 o (x cos <f> 0 + y sin <f> 0 ) + z cos 0 o ] (4 .30) 


E 0 — xE ox 4“ yEoy -f* zE oz 


and 
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Anolytlcal Approximations to the Self-Cell Element 




Figure 4.1: Evaluation of the self-cell element using approximate integration ( — ), 
four-term Taylor series expansion (• • •), and circular disk approximation 
(---)• 
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with 

Eqx = i?o(cos acos 0 O cos <f> Q — sin a sin 4>o) 

Eoy = i?o(cos a cos sin <^ 0 + sinacos<^ 0 ) (4.31) 

Eq z = —Eq cos a sin 0 o 

where a represents the polarization angle of the incident field. It is the angle between 
E‘ and 6. In particular, when a = 0 then H* 2 = 0, corresponding to H-polarization, 
and a = 7r/2 then E\ = 0, corresponding to E-polarization incidence. Upon evalua- 


tion of the current K, the scattered field is given by 

„-jk 0 r 

E*(0, <f>) = -jk 0 Z- N t (0, <f>) (4.32) 

4 itr 

where (r, 6, <j>) are the spherical coordinates of the observation point. Also, 

N t (M) = 0N tg (O,<j>) + (4.33) 

N t $(9, (j>) = cos 6 [cos <f>S x (0, <{> ) -f sin^S y (0, <f > )] - sin 6S X (0, <f > ) (4-34) 

4>) = ~ sin <f>S x (0, <f>) + cos <t>S v (0, (j>) (4.35) 

and 

S(0,<t>) = J jKix^y^^^'^+y'^dx’dy 1 (4.36) 


The field E' can also be described as that attributed to the radiation of the plate 
currents and is responsible for the radar cross section of the plate defined as 

* = ( 4 . 37 ) 

in which p T is a unit vector denoting the polarization of the receiver. 
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First, it is of interest to examine the current distributions on the rectangular plate 
as it has a rather unique and predictable behavior, particularly for principal plane in- 
cidences. Figure 4.3 depicts three-dimensional views of the co- and cross-polarization 
currents on a 2A x 2A conducting plate. An important observation with regard to 
these plots is the high current density values near the edges and the dominance of 
the co-polarized current component relative to the cross-polarized component. The 
singular behavior of the K x currents at the edges is generic to perfectly conducting 
structures with sharp edges. These singularities are responsible for the diffracted 
fields and are the primary source of difficulty in numerical simulations. As 6 in- 
creases, the strength of the cross- polarization currents also increases effecting the 
behavior of the co-polarized currents, especially those toward the back edge of the 
plate. When 0 = 90°, K y have their greatest strength. They are concentrated near 
the side edges and are responsible for the travelling edge waves which, although not 
radiating at backscatter, are crucial in determining the back edge co-polarized cur- 
rent behavior. The lobing structure of the edge currents is particularly interesting 
and unique to all rectangular plates regardless of their size. Generally, for an nA x nA 
plate, the magnitude of the co-polarized currents are associated with n maxima near 
the front and back edges, whereas the cross polarized currents have 2n maxima near 
the side edges. 
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(a) 


Figure 4.3: E-polarization plane wave scattering from a 2A x 2A conducting plate at 
normal incidence (0, = 0°, fa = 90°, a = 90°); No. of samples: 66 x 66; 
FFT pad size: 512 x 512. (a) Co-polarized component of the current 
density, (b) Cross-polarized component of the current density. 
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Since the RCS of a structure is an easily measured quantity, it provides a means 
for validating the solutions. Using the computed plate current densities, the radar 
cross section (RCS) of the plate can be found in accordance with (4.37). Figure 4.4 
illustrates the convergence of the far zone scattered field (using pulse basis) by a 
square perfectly conducting 2A x 2A plate as the size of the FFT pad is progressively 
increased. In all cases the algorithm was terminated when the residual reached a 
normalized value of about 0.01. Also shown in Figure 4.4 is the improved result 
using the CGDFT method. It is observed that an FFT pad of order 1 (minimum 
size) is sufficient when using the CGDFT to yield results that are in agreement with 
the measured data. In contrast, at least a pad of order 3 (along each dimension) is 
required to obtain acceptable results when using the CGFT method and although the 
general behavior of the backscattering cross section approaches that of the measured, 
the convergence to the measured value is not clear near grazing incidence even with 
higher order pads. The principal plane backscatter RCS patterns as computed via 
the CGDFT for the 2A x 2A square plate are compared with the measured data in 
Figure 4.5. The results are seen to be in very good agreement in this case. 

Often of interest is the computation of the plate’s edge-on scattering. As is 
well known, for edge-on incidences the plate currents are rapidly varying and this 
makes their computation a more challenging task. The accuracy of the proposed 
formulations can therefore be best evaluated by examining the edge-on scattered 
field. Some measured data for the edge-on radar cross section have been reported 
in the literature. For example, Figure 4.6 shows the edge-on behavior of a plate 
of constant width ( b = 2A) and varying length (2A < a < 2.5A) with the electric 
field aligned with the shorter side reported in [51]. As can be seen, they compare 
quite favorably with corresponding values computed via the CGDFT formulation. 



Radar Cross Section/A*. dB 
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Angle of incidence 0© (deg.) 


Figure 4.4: Comparison of backscatter RCS patterns for a square 2A x 2A conducting 
plate as computed via the CGFT and CGDFT methods using FFT pads 
of various orders (E-Pol., normal incidence). 
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It should be noted that the CGFT method employing the continuous transform of 
the Green’s function was found inadequate for an accurate prediction of the edge-on 
scattering behavior [48]. 

The radar scattering from a polygonal plate is shown in Figure 4.7 along with 
the corresponding measured data [52]. The scattering characteristics of geometrically 
complex targets may also be simulated by approximating the target by a polygon 
of n sides. This is illustrated in Figure 4.8 where the plate has been modeled as a 
polygon of 8 corners. 






Radar Cross Section a/X*, dB 


90 



Side Length, a/X 

Figure 4.6: E- polarization edge-on RCS for a rectangular metalic plate. 



Figure 4.7: E-polarization scattering from an irregular-edged conducting plate, (a) 
Geometry, (b) Elevation- plane backscattering CGDFT result ( — ) com- 
pared with measured data [52] (o o o). 





Radar Cross Section cr/A*, dB 



ASPECT ANGLE (DEGREES) 

Figure 4.8: E-polarization scattering from an irregular-curved conducting plate, (a) 
Geometry, (b) Near-grazing conical ( 6 = 80°) backscattering CGDFT 
result ( ) compared with the MoM data ( — ). 
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Resistive plates are considered next. In practice conducting surfaces are replaced 
with resistive cards for cross section reduction purposes. By defining the surface 
resistivity Z a as a function of position, different resistivity tapers may be treated. 
As an example, the resistivity can be expressed by a nonlinear function 


Z a (x,y) — Z c + (Z. — Z c ) 


Y Jizfgl V , py-^l V > 

A *f/2 ) \ Vil 2 ) 


(4.38) 


where Z c and Z e may be considered as the resistivities at the center and the edges of 
a rectangular plate, respectively and r x and t v denote the tapers in the corresponding 
directions. Figure 4.9 shows the effect of uniform and non-uniform (parabolic) resis- 
tive tapers on the monostatic cross section of a 2A x 2A plate. The bistatic behavior 
of a polygon of 5 sides in the azimuth plane of 6 = 60° is shown in Figure 4.10. Also 
shown is the result for the same plate with a parabolically tapered resistivity given 

by 


Zo 

(a;/ A) — 2 

2 7 

[(y/A)- 2] 

2 

2 

+ 4 

2 


(4.39) 


where Z 0 is the intrinsic impedance of free space. 

The convergence characteristics of the CGDFT solution for a square 2A x 2A 
conducting plate illuminated under normal incidence is shown in Figure 4.11. At 
each iteration, both the normalized residual error, R t and the incremental error in 
the backscattering cross section, A a are given. These are respectively defined as 


Rs 


M(Jm] ~ E'|| 
Ill'll 


(4.40) 


and 


Aa = a m - , dB (4.41) 

where m denotes the iteration number. It is observed that for far field calculations, 
accurate results can be obtained in much less number of iterations than that required 
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to reach the true solution for the current density. In this case for example, an 
incremental error of 0.1 dB (a relative error of 2%) in the backscattering RCS was 
reached within only 6 iterations. At this point, the normalized residual error was 
about 18%. 



Radar Cross Section a/X*. dB 
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Angle of Incidence 6, deg. 


Figure 4.9: E-polarization scattering from a square 2A x 2A plate with and without 
resistive taper, (a) Backscatter RCS patterns for conducting plate ( — ), 

uniformly tapered plate Z, = Z Q j \ ( ), and parabolically tapered plate 

(■'■)• (b) Three-dimensional view of the parabolic resistive taper for the 
plate (Z. min = 0). 
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(b) 
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Figure 4.10: E-polarization scattering from an irregular conducting plate with and 
without resistive taper. Incidence angles: 0 O = 60°, <j> 0 = 0°; Sampling 
density = 225/ A 2 , (a) Geometry, (b) Bistatic scattering patterns with 
(- - -) and without (— ) resistive taper, (c) Three-dimensional view of 
the parabolic resistive taper for the plate. 
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Scattering from a Polygonal Plate 

Number of Unknowns* 4100, Conical Cut 0*60*. E— Pol. 



(b) 




Normalized Residual Error, R, 
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Number of Iterations 


Figure 4.11: Convergence rate of the CGDFT solution for a square 2A x 2A conduct- 
ing plate at normal incidence. 


Incremental Error in o Q (dB) 
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z 



Figure 4.12: Geometry of an arbitrarily oriented Hertzian dipole in the vicinity of a 
plate. 

4.2.4 Radiation of a Dipole in the Presence of a Plate 


In this section, we consider the problem of radiation by a Hertzian dipole in the 
presence of a resistive plate, illustrated in Figure 4.12. The dipole is centered at 
(sijJ/i, *i)> is of length l <C A and carries a constant excitation current equal to 
unity. Its presence excites currents on the resistive plate which contribute to the 
overall radiation pattern. To compute the plate currents we must solve either (4.12) 
or (4.18) with the incident field given as 


E' x = (Erf + E e d) - x (4.42) 

Ey = ( Eri? + EgiO') • y (4.43) 


where the primes indicate spherical system parameters measured with respect to the 
coordinate system at the dipole as shown in Figure 4.12. We have 

g-ifc.r' 


E r > — 2 Z 0 k 0 t 
E$> = jZ 0 k 0 £ 


4?rJb a r /2 






(4.44) 


(4.45) 
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in which z 1 denotes the dipole orientation and can be represented as 


z 7 = cos 4> r sin e T x + sin <j) T sin <f> r y + cos 0 r z (4.46) 

where ( 0 T ,<f> r ) are the spherical angles of the dipole axis with respect to the plate’s 
coordinate system. Also, r' = r'r 7 is the vector drawn from the dipole’s location to 
the observation point on the plate and 


6 ' = ^ 




x F| 


x f 


(4.47) 


Numerical results based on a solution of (4.12) are given in Figures 4.13 through 
4.18. An FFT pad of order p = 2 and piecewise sinusoidal basis functions were 
used to generate these results. In particular, Figure 4.13 illustrates the current 
components on a perfectly conducting and a resistive 1A x 1A rectangular plate due 
to illumination by a vertical electric dipole positioned at a distance A/ 4 above the 
center of the plate. Figure 4.14 shows the spectrum (the magnitude of the Fourier 
transform) of the x-component of the current density. The principal plane radiation 
patterns of the dipole are shown in Figure 4.15. As seen, the pattern computed with 
the CGFT is in good agreement with that based on the MoM. The corresponding 
results for a horizontally oriented dipole above the same plate are given in Figures 
4.16 through 4.17. 

Finally, the improvement obtained in the convergence rate of the CGFT technique 
when using piecewise sinusoidal (PWS) surface expansion functions is illustrated in 
Figure 4.18. As before, an estimated 100 percent improvement in the convergence 
rate was observed when the PWS basis functions were employed. 




Figure 4.13: The excited surface current density on a 1A x 1A conducting plate ir- 
radiated by a vertical Hertzian dipole positioned a distance A/4 above 
the center of the plate; 25 x 25 unknowns and FFT pad of order = 2. 
(a) X-component. (b) Y-component. 
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Figure 4.14: The spectrum of the surface current density (x-component) on a 1A x 1A 
plate irradiated by a vertical Hertzian dipole positioned a distance A / 4 
above the center of the plate (25 x 25 unknowns and FFT pad of order 
Q — 2 )- 



(a) 


(b) 


Figure 4.15: The principal plane radiation patterns of a Hertzian dipole vertically 
positioned at a distance A/4 above a 1A x 1A conducting plate computed 
by the MoM ( ), and CGFT (• • •) using sinusoidal basis functions, 

(a) Eg pattern, (b) E $ pattern. 
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Figure 4.16: The excited surface current density on a 1A x 1A conducting plate irra- 
diated by a horizontal Hertzian dipole positioned a distance A/4 above 
the center of the plate; 25 x 25 unknowns and FFT pad of order £> = 2. 
(a) Co-polarized component, (b) Cross-polarized component. 



Figure 4.17. The principle plane radiation pattern ( E $ ) of a Hertzian dipole hori- 
zontally positioned at a distance A/4 above a 1 A x 1A conducting plate 
computed by the MoM ( — ), and CGFT (• * #) using sinusoidal basis 
functions. 




Normalized Residual Error, R 
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CGFT Convergence Pattern 



Iteration Number 


Figure 4.18: Improvement in the convergence rate of the CGFT solution for the 
problem of a hertzian dipole positioned above a conducting plate: Si- 
nusoidal basis functions (•••), conventional FFT (delta basis) ( — ); 
25 x 25 unknowns and FFT pad of order g = 2. 



110 



Figure 4.19: Geometry for a dielectric cylinder illuminated by a plane wave. 

4.3 Scattering by Dielectric Cylinders 


We now turn our attention to the problem of scattering by an inhomogeneous 
isotropic dielectric cylinder of relative permittivity e r , as shown in Figure 4.19. The 
cylinder axis coincides with the z-axis and is infinite in extent in this direction. It is 
illuminated by a plane wave given by (4.28) incident at an angle 6 = v/2 

E* = E 0 e^°( x cos + y sin <M (4.48) 


where 


E 0 = £ 0 [sin a(— x sin <f> 0 + y sin <f> 0 ) — z cos a cos <f> 0 } 

To solve for the scattered field, we Introduce the equivalent polarization current 
density (see (3.23)) 

J = - 1)E T , (4.49) 

within the cylinder, where E r is the total field given by 


E T = E* + E 8 = Z C J, 


(4.50) 



Ill 


in which 


2c = 77 


(4.51) 


jk 0 {e T - 1) 

The governing system of integral equations is now obtained by substituting for the 
scattered field in (4.50). For an arbitrary polarization of incidence we have 

E \p) = Z c (p)J(p) + jk 0 Z 0 jf J(p') • T( P] p')ds (4.52) 


where T is now given by 

/ , 1 d 2 . 

L2 ^ 


r = 


i d 2 


k 2 <9z 2 ’ k 2 dxdy 


i a 2 

Ar 2 dydx 


1 S 2 


G c (p; /»') 


(4.53) 


and 


4j 


(4.54) 


Following the same procedure discussed earlier in connection with the plate prob- 
lem, we may write 


M-lAT-l 

= } ] } ^ Jmn * !/) 

m=0 n=0 


(4.55) 


where 


$mn(x,y) = - * m ,y - yn) ( 4 - 56 ) 

$(z, y) = zz<£ x (z, y) + yy<j> v (x , y) + zz<f> z (x, y) (4.57) 

and 4> t denotes the expansion function in the ^-direction. 

The two-dimensional continuous Fourier transform of G c is given by (Appendix 

B) 


1 


G c (k x , k y ) — £ 2 ^ 


(4.58) 
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However, as noted earlier, the use of (4.58) in the CGFFT solution of (4.52) will 
result in aliasing errors. An additional difficulty will also arise because of the ring 
singularity of (4.58) occuring when kl + k* = k\. The inherent approximation in 
the implementation of the inverse FFT that the transform be constant over each cell 
is obviously not valid for those cells coinciding with the ring singularity. This can 
cause substantial errors and often leads to the failure of the discrete system as an 
acceptable representation of the continuous one. 

To correct both of the above sources of error, the procedure described earlier 
can be employed here as well. That is, the original continuous integrals are first 
discretized before proceeding with their evaluation via the discrete Fourier transform. 

The new discrete system of equations is expressed as 


E i , J = Z ctJ J 0 + ^^A,,-J r 


(4.59) 


where 



C xv 

£ 1 * 



\ 


0 


o C x ) 


(4.60) 


which should be compared with (4.15) and (4.16). Similarly, the non-trivial entries 


of A are now given by 


(if = Hh 2) (x,y; x\y')4> x {x' - x m ,y' - y n )ds' 

(if = #o a) (*»v;*'»v')^(*' ~ x m ,y' - y n )ds' 

(if = H^\x,y\x\y')<i> x {x' - x m ,y' - y n )ds' (4.61) 

(if = (1 + jj Ho\x,y\x',y')4> v {x' - x m ,y' - y n )ds' 

(if = ! 4 2) G(x, y; x', y')4> t {x' - * m , y' - y n )ds' 
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which may be considered as the ‘cylindrical counterparts’ to (4.17) since <r mn here 
denotes the incremental surface element corresponding to the mnth cylindrical cell 
on the target. Again, these expressions are evaluated at (x,y) = (x,, yj). Applying 
the convolution theorem to (4.59) yields the final form of the system of equations 


where 


If 7 ~ 

E,j = ZdjJij + DFT -1 {A • J} 


A H 


kl-Dl -D X D V 0 
-D y D x kl-Dl 0 


(4.62) 


(4.63) 


and D x and D v are given by (4.21)-(4.22). In the above, C is the discrete Fourier 
transform of the sequence 


Cm. = / (4.64) 

J0mn 

This integral can be evaluated numerically except when (m = n) which corresponds 
to the self-cell term. When (m = n) we must resort to an analytical evaluation 
similar to (4.26) as [36]: 

Coo a J i; [irir 0 H 1 <1) (tr„)-2j) (4.65) 


with r 0 as given in (4.27). 

To illustrate the accuracy of the above formulation, the bistatic scattering from an 
infinitely long triangular cylinder is shown in Figure 4.20. The cylinder is perfectly 
conducting and as seen, our result agrees very well with a corresponding direct 
solution provided in [53]. 



0.7071 X 



20 1 

00 



Angle of Incidence 0, deg. 


Figure 4.20: Scattering from a conducting triangular cylinder illuminated by an E- 
polarized plane wave, (a) Geometry, (b) Comparison of bistatic echo 
widths obtained from the CGDFT method ( — ) and the direct method 
(• • •) [53]. 
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4.4 Summary 

Two formulations for a conjugate gradient solution of the scattering by plates 
of arbitrary shapes were presented. One of the formulations (CGFT) employed the 
sampled (truncated) continuous transform of the Green’s function for the evaluation 
of the convolution integrals. The other (CGDFT) employed finite duration discrete 
Fourier transforms for the evaluation of the same integrals. As with the strip problem 
studied in the previous chapter, the latter method was found to provide an accurate 
simulation of the plate scattering by eliminating aliasing errors (other than those 
due to under-sampling). It was also found to be substantially more efficient than the 
former method. Furthermore, it was noted that the convergence of the solution is 
substantially faster for plates of non-zero resistivity and this is attributed to the less 
singular behavior of the edge currents. 

The CGDFT method was also applied to the problem of scattering by dielectric 
and conducting cylinders of arbitrary cross sections and a degree of accuracy and 
efficiency similar to the plate problem was observed. 



CHAPTER V 


GENERALIZED IMPEDANCE BOUNDARY 

CONDITIONS 


5.1 Introduction 

Generalized Impedance Boundary Conditions (GIBC) are higher order bound- 
ary conditions which involve derivatives of the fields beyond the first. They have 
been found to be more effective than the traditional first order (standard) conditions 
(SIBC) in modeling thick dielectric coatings and layers. The GIBCs offer several ad- 
vantages in both asymptotic and numerical analysis of electromagnetic problems. For 
example, in the case of high frequency analysis, they allow an accurate replacement 
of a coating on a layer with a sheet boundary condition amenable to a Wiener-Hopf 
analysis [14, 54], or some function theoretic approach [55]. In numerical analysis, the 
profile of a coating can be replaced by a simple boundary condition on the surface of 
the coating. This eliminates the need for introducing unknown polarization currents 
within the volume of the coating or material layer, thus leading to a more efficient 
solution from the numerical standpoint. 

A convenient form of these conditions is expressed in terms of the normal deriva- 
tives of the field components 
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M 


flm d m E n 

bo dnm 

o' STH n 


M 


= 0 


= 0 


(5.1) 


(-J k o) m dn' 

where the subscript n denotes the normal conipoiient to the surface and o m and ci m 
are constants specific to the material and geometrical properties of the structure. 
These constants are chosen so that the application of the boundary conditions will 
reproduce the desired scattering behavior of the surface or coating layer under con- 
sideration. They can be derived by employing a suitable expansion of the coating’s 
Fresnel reflection coefficient and by matching the constants of the expansion to those 
of the compatible conditions implied by the GIBC (5.1). Since for a given problem, 
the electric and magnetic fields may not be specified independently, the coefficients 
a m and a' m are related by the relation [13] 


J2 Om 

l m=0,2,— 


E «:) = ( E •-)( I 

n=0 t 2,»* / \m= 1,3,-* / \m= 1,3, * 


(5.2) 


which is a form of duality condition. The integer M in (5.1) is the order of the 
conditions. For example, when M = M' = 1, we have a first order condition 


£„ - — ^ = 0 


H. - 4-^=0 


(5.3) 


jk 0 rj dn ’ " jk 0 dn 

The above first order condition can be shown to be equivalent to the standard 
impedance boundary condition (SIBC), also known as the Leontovich boundary con- 
dition [56]- [58], provided r? = a 0 /a x = a[/a' 0 is identifier as the normalized surface 
impedance of the sheet. The SIBC is often used for simulating material coating on 
conducting bodies and in this case the normalized surface impedance tj is given by 


rj = — = j— tan (Rk 0 t) 

a\ Cr 


(5.4) 
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where t is the coating thickness and N = y/e r fi r is the index of refraction. The validity 
of this SIBC has been examined by several authors [57], [59]-[61] and in general it 
holds when 

|K| > 1 (5.5) 

and 

|SmN|M > 1 (5.6) 

The first condition ensures that within the medium, the field behaves essentially 
as a plane wave propagating in the direction of the inward normal to the coating. 
The second condition, on the other hand, imposes the requirement that the inward 
traveling field suffers enough attenuation so that no outward traveling waves exit at 
the interface due to reflection. Abo, for inhomogeneous materiab, the SIBC remains 
valid if the lateral variations of the impedance in the medium are slow, that is 

« 1 < 5 - 7 > 

where V denotes gradient in coordinates transverse to the normal. 

Inherently, the SIBC does not permit modeling of polarization currents which are 
normal to the layer and is thus most suited for near normal incidences. However, by 
increasing the order of the condition, it is possible to allow accurate simulations of 
fairly thick layers and unlike the SIBC, the accuracy of the higher order conditions 
improves as the incidence angle approaches grazing. 

A third order GIBC derived recently [62] for the simulation of high contrast 
dielectric coatings is given by (5.1) with M = 3 and 

- = (*- a) -*»(£)] 
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Oi = ~Jt r 


1 + tan(KA: < ,<) tan 


(£)] 


and 


a 2 


= If (K) + k °‘ ( N - |f) 

1 + tan(Nfc<>0 tan j | 


<*3 — 


j k 0 ttr 

2N 


tan(Nfco<) 


" tan (w)l 


< = ( H -^) 


1 + tan(NA: 0 t) tan 


@ 9 ] 


ai = jftr tan(NM) - tan 


(5.8) 


= ^ { 1 + tan (KkJ) tan ) - k a t (« - -^) (5.9) 

• jtan(N* 0 i) - tan j j 


a-. = 


jkptfl r 

2N 


1 + tan(tt& 0 t) tan 


(*)] 


The GIBCs are usually applied at the top layer of the surface under study. How- 
ever, in some applications it is desirable to apply the conditions at another plane 
of reference. This is convenient for a coating on a ground plane where one may be 
interested in invoking the image theory. In such cases the original coefficients a m 
must be replaced by (Appendix C) 

= f; « m -n . m = (5.10) 

n= 0 n * 
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where t denotes the transfer distance. 

5.2 Two-dimensional Impedance Inserts 


In the two-dimensional case, the impedance insert is assumed to be infinite in 

Q 

length with no field variations along the z-direction (— = 0). The insert is assumed 
to satisfy a generalized impedance boundary condition at its surface. This may serve 
as a model of a partially coated conducting plane or a material filled groove whose 
scattering behavior under plane wave illumination is of interest (Figure 5.1). 
Introducing the equivalent magnetic current density over the insert we have 


M = E x n 


M z = E x 
M x = -E z 


(5.11) 


and by imposing the continuity of the tangential field components an integral equa- 
tion for the current density can be obtained and solved numerically. 

In the following, we will derive integral equations for the equivalent current den- 
sity based on a third order generalized boundary condition (M = 3). The two 
principal polarizations axe treated separately. 


5.2.1 H-Polarization 

The incident fields are assumed to be of the form 


H’ 

_ £ e «* 4>o+V * 0 ) 

(5.12) 

E* 

= Z 0 (x sin 4> 0 -y cos <j> 0 ) e jk ° {xco •*.+»*»«•> 

(5.13) 


and the first of the conditions (5.1) is relevant in this case. For a third order condition 


f °m STE, 
£lt (-A)" 


we write 
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X 

Figure 5.1: Simulation of a partially coated conducting ground plane by an 
impedance insert. 
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or equivalently, 

( 0.2 d 2 \ ( a i a 3 d 2 \ d 

{'-swn+inz+iswh *- 0 (5l4) 

In order to derive an integral equation on the basis of (5.14), it is desirable to 
work with transverse derivatives. This allows for a convenient application of the 
Fourier transform to solve the resulting integral equation. To do so, we note from 
the divergence relation that 


V-E = 0, 

and from the wave equation 


3Ey 

dy 


dE x 

dx 




-(*•■*£) 


E„ 


Introducing these into (5.14) along with (5.11), we have 

1 d 2 


1 4- E 2 1 1 + 


k 2 dx \ 


(5.15) 


(5.16) 


, £ ' + ^r + ^( 1 + w)]|v.=° (5.17) 


where 


F t = — , £=1,2,3 (5.18) 

do 

and E y is the component of the total field normal to the coating. It can be expressed 
as the sum of the geometrical optics field in the absence of the sheet (short-circuited) 
and the scattered field in its presence 

E y = E?° + E’ y = Ei + E r y + El (5.19) 

Since the tangential electric field vanishes over the conducting ground plane, the 
reflected field is given by 


E[ = E' y = -Z 0 oos<f> 0 ei k ° xcot +°, 


y = 0 


(5.20) 
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Also, the scattered field can be expressed as as 

d 


Ey — “T [ W/ 2 2M z {x')-2-H? ) {k 0 \x - x'\)dx' 
v 4 J J—w/2 ox 


( 5 . 21 ) 


where the factor of two has been introduced in accordance with the image theory. 
Substituting these into (5.17) yields 


/ 1 d 2 \ 

r 

1 +F2 V + k% dx 2 ) 

\ 4 j dx J- 


(5.22) 


= 0 


- Z. cos ***--*•} + Jg [* + A (> + . 

To eliminate the x derivative, we integrate both sides with respect to x and obtain 

Fi + F 3 (l + Jij-i)] M >< x ) 


(l + Fj sin 3 fa) Z 0 e? k ° x<aa *° = - 


(5.23) 


+ 


1 + 


Fi (‘ + i* J?)] J-/2 x '\ )dx ' 


which is the desired integral equation in M z . 


5.2.2 E-Polarization 

The incident fields in this case are of the form 


jko|iOM^+»*i«*o) (5-24) 

H* = -(xsm<f> 0 -ysin<f> 0 )Y o e jk ^* o+v ^ o) (5.25) 

and the appropriate boundary condition is given by the second relation in (5.1) with 
M = 3 




+ A£.) 



(5.26) 
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Again, using the divergence relation and the wave equation and following steps sim- 
ilar to those taken for H-polarization, we obtain 


^ +r ’( 1+ iJ9 


H x + 


’HS1 


I 1 + F 2 ( 1 + isihs ) I y oM - = 0 (5.27) 

in which the definition of F' t is analogous to that given for Ft (see (5.18)) and 

Hr = H°° + H’ = H' X + HI + HZ (5.28) 

Imposing the condition on the tangential electric field, we find that 


H r x = H x = —Y 0 sin <f> 0 e jkoXCO **‘° , y = 0 


(5.29) 


and write the scattered magnetic field as 

= Cr (‘ + hl&) / - AW ■ (5.30) 

Substituting these back into (5.27) yields the integral equation 

(fj -fusin’ A,) = i M,(x) 


F,+H 1 + 


1 d 2 ' 
kj. dx 2 


7 t 1 + w) Cl - x'Drfi' 


(5.31) 


to be solved for M x . 


5.2.3 Specialization to SIBC 

As mentioned earlier in this chapter, the SIBC formulation has been traditionally 
applied to coated conducting bodies as well as dielectric filled metal-backed cavities. 
The above integral equations can be readily reduced to those corresponding to the 
SIBC formulation by setting Ft = F' t = 0, l = 2,3 in accordance with conditions 
(5.1). Doing so yields 

+ J Cl M,(x')HW(K\ x - x'l W 


(5.32) 



125 


for H-polarization and 


sin + 7 (l + 1 | J^) Cl ~ 


(5.33) 


for E-polarization where ^ = F x = l/i^' = ai/a c is identified as the normalized 
surface admittance of the insert. 


5.2.4 CGFFT Formulation 

The integral equations (5.23) and (5.31) are amenable to a solution via the 
CGFFT method. To put them into a suitable form, we first discretize the magnetic 
current density using the piecewise constant basis functions and follow an analysis 
similar to that presented in Section 3.3.4 in connection with the scattering from a 
strip. The integral equations axe then put into the form 

[l + F 2 (x) sin 2 <j> 0 ] Z 0 e jk ° xco '*° = 

\ F,(x)M z (x) + ±F 3 (x) DFT - 1 {( k\ - D 2 )M,}] 

+ DFT - 1 {M z f } + ^F 2 (x) DFT - 1 {(* 2 - D\)M Z f } (5.34) 

and 

>,(*) + r 3 (x) sin 2 *.] sin 4 >0 e‘ k - r ~*’ = 

i M.(x) + ifi(x) DFT {(* 2 - D 2 )®,}] 

+ EM dft -■ {(fc* _ Dl)M,r} + ^ DFT -* {(A : 2 - } 


(5.35) 
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where the spatial dependences of the F and F' coefficients are expressed explicitly 
to allow for slow lateral variations in the electrical properties of the sheet and T 
is the discrete transform of the train T on defined in (3.85). It is noted that the 
spatial derivatives are carried out relatively easily in the transform domain as was 
the case in the strip problem. A CGFFT implementation of (5.34) and (5.35) is a 
straightforward task. 

5.3 Three-Dimensional Impedance Inserts 

In the previous section, we presented an implementation of a third order GIBC 
for scattering by a two-dimensional impedance insert simulating an infinite groove 
recessed in and/or a coating on a ground plane. Here we present a corresponding 
implementation for the three-dimensional case. 

5.3.1 The Integral Equations 

Consider the geometry shown in Fig. 2 illuminated by a harmonic plane wave 

H‘ = H 0 e"-' M ** r) (5.36) 

E‘ = Z 0 H’ x k { (5.37) 

where fc, is given by ( 4.30) 

ki = - [sin 6 0 (x cos <p 0 + y sin <f> 0 ) + z cos 6 a ] 

k 0 is the free space wave number, and Z 0 = l/Y 0 is the free space intrinsic impedance. 
Also, 

H ox = y o (sin a cos 6 0 cos <f> 0 -f cos a sin <f > 0 ) 

H oy = y[>(sin a cos 9 0 sin <f> 0 — cos a cos <f> 0 ) (5.38) 

H oz = —Y 0 sin a sin 6 0 



127 


and 


E ox = cos a cos B 0 cos <f> 0 — sin a sin <j> 0 

Eoy = cos a cos 6 0 sin <f> 0 + sin a cos <f> 0 (5.39) 

E ox = — cos a sin 6 0 


in which a represents the polarization angle of the incident field (when a = 0 then 
Hi = 0, corresponding to H-polarization, and when a = tt/ 2 then E' z = 0, corre- 
sponding to E-polarization incidence). 

As before, the application of the conditions (5.1) over the surface of the impedance 
sheet requires the introduction of a magnetic current density vector M(x, y) as de- 
fined in (5.11) with both transverse components present. A surface integral equation 
for M can then be derived by following a procedure similar to that discussed for the 
two-dimensional analysis. Before doing so, however, it is again instructive to rewrite 
the boundary conditions in terms of the tangential derivation. This is expected to 
directly yield a symmetric set of equations with respect to M x and M v . From [13] 


we find that conditions ( 5.1) are equivalent to (M = 2) 


- PZH + -Lo^ + —^z^ 

- —PZ 0 H V + jk Q dx + jko p,A> Qy 


(5.40) 


E y 


PZ 0 H X + tt" 
3 * o 



jk 0 P' ° dx 


where 


P 

Q 


02 + Op 
O l 

a 2 
a i 


(5.41) 

(5.42) 


and analogous expressions hold for the primed quantities P* and Q . To derive the 
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integral equation, we invoke (5.11) (note the new coordinate system) 


M = E x n 


and rewrite the conditions (5.40) as 


M X = Ey 


My = -E X 


pM, Z,H„+ , k ' P di + ]ko pp, z ° g y 


J _ _L_QL 7 9h ?° 

° 11 jk,P dx jk.PP' ° By 


— M — Z H‘ 1 Q , \ Q ^ 

P 1 ° * jk 0 P dy jk c PP' ° dx 


= 7ff ao , l Q dE?° 1 Q „ dH?° 

° * jk a P dy jk 0 PP' ° dx 


where we have also made use of the definitions 


E = E* + E r + E* = E ao + E* 


(5.43) 


(5.44) 


H = H‘ + H r + H' = H ao + H' 


Again, the superscript GO specifies the geometrical optics fields in the absence of the 
sheet and the superscript s specifies the scattered field in its presence. By imposing 
the boundary conditions on the tangential components of the electric field over the 
perfect conductor, we find that 


E r = (~E OI x - E^y + E oz z)e~ lk ^ 


(5.45) 


H r = (H 0X x + H^y - H ot z)e~ jk ^ 
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where 


k r = - [sin 0 o (x cos <f> 0 + y sin <f> 0 ) - z cos 6 a ] 


(5.46) 


is the unit propagation vector for the reflected fields. The scattered magnetic field 
may be expressed in terms of the equivalent magnetic current density as 

H' = -jk 0 Y o JJ T(n r') • [2M(r'))^' (5.47) 

where 5 denotes integration over the surface of the sheet and 

F= (l+^Vv)G p (r;r') (5.48) 

is the free space dyadic Green’s function with the factor of two accounting for the 
presence of the ground plane. More explicitly, 

d 2 


h -=- 2 J kL H*'- ^ (*» + J?) + y,) 


dxdy\ 


G 0 {r ; r')ds' 


(5.49) 


H‘v = 


2jY 0 


L 




(« + £)] 


G„(r; r')ds' 


(5.50) 


and 


HI = j $ [M r (*W)^i +M v (x',j/')- — 


Golds' 


'dxdz • " 'dydz\ 

Rewriting the last equation and making use of the identity (Appendix D) 

A G (r;r') = -i«( r-r') 
it follows (from distribution theory) that 


ti» ijf 

Hz ~~r 

H'O 


£m,(x, V ) + 


(5.51) 


(5.52) 


(5.53) 
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To formulate an integral equation from (5.43), we also need to express the normal 
component of the scattered electric field E* z in terms of M. We have 


E-. = lf s S ')|j - M,(x\ ,') 


a_ 

dx 


Golds' 


Substituting now for the field quantities in (5.43) yields 


-J-m x 

2 P 


(5.54) 


+ 


K Is (*• + dx*) + 


d 2 


dxdy 


G 0 ( r; r')ds' 


f 

jk 0 P Js 




^G 0 (r;r')ds' 


+ 


1 Q' d 

2klPP'dx 


d d . . 

— x M,{x,y) + Ty M,(x,y) 


= ^ZoH ox + — sin 0„ sin exp{jk 0 [sin 0 O (x cos rj> 0 + ysin ^ 0 )]} 


(5.55) 


— M y (x,y) 


+ 


t Is + (*’ + $) ] G »< ri ^ 


+ 


J_v f 

jk 0 P Js 


m x (x ,y )^ ^ v ^ x ' y ^dx dx Go{T ' r ' )ds ' 


+ 


i <?' a 

2klPP'dy 


d d 

—M x (x,y) + —M v (x,y) 
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ZoHoy - J sin e 0 cos <f> 0 E 0z ] exp {jk a [sin 0 O (x cos <f> 0 + y sin <f> 0 )] } 


(5.56) 


which is a coupled pair of integral equations to be solved for M x and M y . 


5.3.2 Specialization to SIBC 

The first order (standard) impedance boundary condition is obtained from (5.41)— 
(5.42) by letting a 2 = a' 2 -+ 0 eliminating the corresponding terms in (5.40) 


E x = -PZ 0 H V 


(5.57) 


Ey — P Z 0 H X 

where P = a 0 ja x - y, may be regarded as the normalized surface impedance of the 
sheet. Expressing these conditions in vector form, we may write 

E — (n • E) n = y,Z 0 n x H (5.58) 


In view of (5.57), equations (5.55) and (5.56) reduce to the simpler pair of integral 
equations 


M r (x, y) + 

2y, 


i J s ["•<*•»> (*•’ + &) + 


G 0 (r; r')ds' = Z 0 H' X 


(5.59) 


Ws My(X ' y)+ k 


U[ M ■<*'»' 


)^r. + W* to (*J + £;) ] ' G.(r, r')ds' = Z„Hi 


(P_ 

dxdy 


(5.60) 


As expected when y t — * oo, indicating an open gap in a conducting screen, these 
equations further reduce to the dual of those pertaining to a perfectly conducting 
plate already discussed in Chapter 4. 
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5.3.3 CGFFT Formulation 

In order to put the integral equations (5.55) and (5.56) into a form suitable for a 
CGFFT, we expand M(x, y) in piecewise constant basis functions. In particular, we 
set 

M(x,y)= £ „P{x - x p )P{y - y q ) (5.61) 

p=0 9=0 

where x v = pAx + , y, = qAy + ^ and M m = xM xpq + yM v „ represent the 

unknown coefficients of the expansion function. Employing this expansion, we may 
rewrite the surface integrals as 

f M(r') • r o (r; r')ds' = ^ (5.62) 

JS p=0 q = 0 

where Z<j is given by ( 4.16). Applying the discrete Fourier transform, the above 
further reduces to 

/ M(r') • r o (r; r ')*' = DFT "'{I • K) (5.63) 

V S 

where S denotes the discrete Fourier transform of S and is given by ( 3.109). Em- 
ploying this result in (5.55) and (5.56), yields 

F 1 (x,y)M x (x,y) 


+ ^DFT- 1 {[M,(l: 0 2 - Dl) - M„D,dA A 

"'O 

+ t|~F 2 (x, y)DFT _1 { (M t D\ - M y D x D y ) £} 


- |f DFT" 1 {M X D\ + M y D x D y ) 
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= 2 [ Z 0 H 0X + F 2 (x,y) sin 9 0 sin <f> 0 E 0Z ] exp {j k a [sin 6 0 (x cos 4> 0 + y sin <f > 0 )] } 

(5.64) 


and 


F 1 (x,y)M y (x,y) 


+ ^DFT- 1 { \-M x D x D y + M v {kl - D\)] |} 
*0 


+ -%rF 2 (x,y)DFT- 1 {(-M x D x D y + M v Dl)t} 
J Ko 


- |fDFT -'{M x D x D y + M y Dl\ 


= 2 [Z o H 0y - F 2 {x, y) sin 9 0 cos <f> 0 E 0X ] exp {jk a [sin 0 O (x cos 4> a + y sin <£„)]} 

(5.65) 


which are applicable for a solution via the CGDFT method. In the above, 



F 2 = 


9 . 

p 


and the spatial dependence of the F coefficients are again shown to indicate the 
presence of a slow planar variation in the electrical properties of the sheet. 


5.4 Summary 


In this chapter, we introduced the generalized impedance boundary conditions 
and studied their incorporation in the general CGFFT formulation. The formulation 
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was applied to the simulation of what could be referred to as generalized impedance 
inserts. In particular, a third order GIBC was applied to the simulation of two- 
dimensional impedance inserts while a second order condition was considered in 
the three-dimensional case. In both cases, the first order (SIBC) formulation was 
obtained by setting the appropriate higher order coefficients to zero. 

The combined GIBC/CGFFT formulations discussed here can be used in the 
study of partially coated conducting planes as well as cavity-backed apertures re- 
cessed in a ground plane. These structures may be adequately represented by 
impedance inserts with appropriately chosen coefficients. This will be the subject of 
Part Two of the Thesis. 
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Part II 

SCATTERING BY CAVITY 
STRUCTURES 


CHAPTER VI 


SCATTERING BY MATERIAL FILLED 

GROOVES 


6.1 Introduction 

The study of electromagnetic scattering from filled cavities recessed in ground 
planes is important in modeling the radar response to various man-made structures. 

In this chapter, an exact full- wave formulation is first developed for the rectangu- 
lar groove problem based on the Generalized Network Theory [63]. This theory has 
been applied to a number of aperture and slot problems in the past [64, 65]. In this 
method the external fields are expressed in terms of the scattering integral while the 
fields internal to the dielectric medium are given in terms of appropriate waveguide 
modes specific to the particular problem. An integral equation is then set up by 
employing the equivalence principle and enforcing continuity of the electromagnetic 
fields across the interface. This method, although rigorous, is computationally in- 
tensive and is limited in application to structures whose electrical size is relatively 
small. Moreover, due to the nature of the formulation, a solution is possible only 
for canonical geometries for which the orthogonal wave functions associated with the 
cavity can be found. 

Next, the formulation is specialized to narrow grooves. Analytical expressions for 
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the equivalent magnetic current distribution over the aperture of narrow grooves are 
derived based on a quasi-static approximation of the pertinent integral equations. 
The solutions exhibit the expected edge behavior at the terminations and are used 
to find closed form expressions for the echo width of the groove. 

Finally, an approximate formulation based on the GIBCs is presented and shown 
to be amenable to the CGFFT method of solution having an order O(N) memory 
requirement. In contrast, the exact integral equation does not lend itself to such a so- 
lution and must be solved by a matrix inversion approach having an (D{N ) memory 
requirement. It is found, unfortunately, that the GIBC formulation yields satisfac- 
tory results only when the contributions of the groove’s terminations are negligible. 
This is because the GIBCs were derived for a coating without terminations and must 
be supplemented by more accurate conditions in the vicinity of such material discon- 
tinuities. A hybrid procedure is, therefore, introduced that combines the exact and 
GIBC formulations. The proposed procedure utilizes the solution obtained from the 
GIBC /CGFFT in a region sufficiently away from the terminations and then finds the 
near-edge currents based on the exact formulation. Despite an increase in the com- 
plexity of the formulation, the memory requirement of the hybrid method remains 
essentially of the order 0{N) and can be used when the material constituency of the 
filling does not allow the application of neither SIBC nor GIBCs of higher order. 

6.2 Full- Wave Formulation 

Consider the infinitely long groove of width to, and depth d illuminated by the 
plane wave 

H‘( or E’) = zeJMsWo+vwn**) 


( 6 . 1 ) 
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Figure 6.1: Geometry of a filled rectangular groove. 

for H- (or E-) polarization, where k a = 2ir/A is the free space wave number and 4> 0 
is the angle of incidence (Figure 6.1). The groove is assumed to be filled with a 
material of index of refraction N = y/t T ^ T . A standard approach to formulate the 
scattered field by the groove is to employ the equivalence principle [66]. Accordingly, 
the aperture is closed by a perfect conductor and the equivalent magnetic current 
(Figure 6.2) 


M = Exn = Exy (6.2) 

is introduced over the aperture at y = 0 + . The scattered fields outside the cavity are 
those radiated by the equivalent magnetic current and consistent with the continuity 
of the tangential electric field, the field inside the cavity is that radiated by -M 
placed at y = 0 across the aperture. To find the equivalent magnetic current we 
must also enforce the continuity of the tangential magnetic field across the aperture. 
We have 


n x [H“(M) + H'“] = n x H*(M') 


(6.3) 
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where H‘° is the total field on the ground plane in the absence of the groove (aperture 
short-circuited), H° represents the tangential scattered field above the aperture and 
H fc is the total field below the aperture. To construct an integral equation in M, 
H° and H fc must be expressed in terms of the Green’s function corresponding to 
each region. The external scattered field (attributed to M) can be expressed as the 
surface integral 

H°(r) = -jk 0 Y 0 r 1 2M(x') • r(x;x')dx' (6.4) 

J- w/2 

where Y a is the intrinsic admittance of the free space and f is the two dimensional 
dyadic Green’s function 

f(x; *') = -j (1 + + H H?\k,\x - x'|) (6.5) 

and a factor of two was introduced in (6.4) based on the image theory to account for 
the presence of the ground plane. 

The internal fields (those attributed to M') can be written in terms of the TM z and 
TEzwaveguide modes as 

E 6 = E™ + E te = -jk b Z b {zV™) - V x {zV TE ) (6.6) 

H 6 = H™ + H T£; = V x (zV™)-jk b Y b (z* TE ) (6.7) 

where fcj, = \$k 0 is the wave number inside the cavity and Zf, = l/Tj, is the intrin- 
sic impedance of the filling material. The functions ty™ and are the wave 

potentials both satisfying the scalar wave equation 

(68) 

subject to the boundary conditions 

E x = E z = 0 
E v = 


E z = 0 


y = —d 
x = ±w/2 


(6.9) 

( 6 . 10 ) 
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on the cavity walls, and 

M' = Exn y = 0 (6.11) 

over the aperture. 

Below, we consider the two principal polarizations separately. 


6.2.1 H-polarization 

For H-polarization (TEzcase) we have 


H ia _ H mc + H T* _ 2e ,fcoX cc*. 


( 6 . 12 ) 


which is the geometrical optics field in the absence of the groove. The tangential 
component of the external scattered field is given by 

jT* M,{x')H ? - x'\ )dx' (6.13) 

while the internal fields are given by (6.6) and (6.7) and in this case we have 

E‘ = E™ = -V x (z*™) = E + 9^ TE (6.14) 

H‘ = H TE = -jkM~z<t TE ) (6.1.5) 

In order to find useful expressions for E 6 and H 6 , we need to solve for the wave 

potential '& TE . To this end, ty TE can be expressed as an infinite sum of orthogonal 
modes 

V TE = A p t/;$ (6.16) 

p=0 

where ipp are the waveguide modes all satisfying the wave equation (6.8) and A p are 
coefficients to be found. Substituting for V TE in (6.14), and using (6.9)-(6.10), the 




H' = z Y gjyx cos<) 0 +y sin(j) 0 ) 


Equivalence Principle: M=Exn 



Figure 6.2: Application of the equivalence principle to aperture problems. 
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boundary conditions to be imposed on the cavity walls are 


irti = o 


dy^ 

d 

dx' rp 


= -d 


o 

— =0 x = ±w/2 


(6.17) 

(6.18) 


and a set of eigenfunctions which satisfy these and the wave equation is 


= cos[Ar p (y + d)] cos [~(x - tw/2)] 


where k p satisfies the separation parameter equation 


k\ = kl - (^)> 


(6.19) 


( 6 . 20 ) 


We now seek to find the coefficients Ap. Upon enforcing the condition (6.11) on the 
aperture ( y = 0), we find that 


T. k p A p sin (kpd) cos [~(* — «»/2)] = —M' z = M z (6.21) 

Multiplying both sides by cos [— (x — tv/2)] and integrating over the aperture yields 

E kp-Ap sin (kpd) f“ cos [— (x — tp/2)] cos [— (x — u;/2)]dx 
„ J-w/2 W W 


J-w/2 


M t (x) cos [— (x — w/2))dx 

-w/2 W 


( 6 . 22 ) 


and by invoking the orthogonality relation 


/ 


cos [— (x — w/2)] cos [— (x — to/2)]dx 

— w/2 W 


W 


W 


■^■[l + f>po] p = q 
3 p ± q 


( 6 . 23 ) 


we find 


A — I 

[l f tvkp sm ( kpd) J — w/2 


f“ M x (x) cos [— (x — w/2)]dx 
J—w/2 W 


(6.24) 
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where 6„ is the Kronecker delta defined in (3.117). Thus, the function 'H TE in 
(6.16) is completely defined and the magnetic field in the internal region may now 
be expressed explicitly as 

o° o pir 

*<*•»> - - jhYt g [1 + S^k, sin(M ) + 1,1 008 l V <: C ~ w/2)] 


f * M,(x') cos [— {x' — ui/2)]dx' 
J-w /2 tV 


which at y = 0 gives the tangential field just below the aperture 

*(*■ o) - - jhYi £ [i + uk t ^(M) 003 if (i - - ,/2)l 


(6.25) 


/ / M x (x') cos [— (s' — u>/2)]dx' 
J-w /2 tv 


(6.26) 


This is equivalent to expressing the tangential internal fields as 


H\{-M x ) = -jk b Y h jH^ M x (x')G h (x; x’)dx’ 


(6.27) 


where the Green’s function is given by (see (6.26)) 


G h (x, x 1 ) = £ r l jr-j . cos [^(x - «,/2)] cos p(x' - u,/2)] (6.28) 

±o £ p wk p tan(A^d) tu u> 


£p = 1 + 6p 


(6.29) 


Substituting (6.12), (6.13), and (6.27) into (6.3) we obtain 

Z 0 e? koZCm *° = 1^ M x {x')H?\k 0 \x - x'\)dx 

4 J-w/2 


^ cos [^(x — iw/2)] r°! 2 , .px , 

+ V La Li L_ii / M,(x')cos[— (x' - u>/2)]dx' 

^ e p WT) hp J-w /2 w 
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where rjh p are the normalized H-mode impedances of the cavity given by 


Vhp = j-r*b tanffcpd) 


(6.31) 


and Zf, = yj n r l t, is the normalized intrinsic impedance of the internal region. Equa- 
tion (6.30) is an exact integral equation to be solved for M z (x). 


6.2.2 E-polarization 

For E-polarization (TM z case), we have 

H' x a = -2 Y 0 sin <j> 0 e ik (6.32) 

and the corresponding tangential scattered fields are given by 
k Y f l fp \ /W2 

= -y (l + ) /_ W2 - x-l )dz' (6.33) 

and 

H x (-M x ) = f M x (x')G e (x; x')dx' (6.34) 

To find the cavity Green’s function G e we note that 

E 6 = E™ = -jk b Zb{z*™) (6.35) 

H 6 = H™ = = (6.36) 

ay ox ' 

Following steps similar to those taken for the H-polarization case, the wave po- 
tential ty™ is expressed as 

*™ = f ]B P 4>; (6.37) 

p=o 

and the boimdary conditions to be satisfied on the walls are 

i’l = 0 y = -d (6.38) 

= o x = ±w/2 (6.39) 
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which are satisfied by chhosing ip p as 


= sin[fc p (y + d)] sin [^(z ~ W2)] 


(6.40) 


and k p is defined in (6.20). Enforcing now the boundary condition (6.11) on the 
aperture, we have 

^k p B pS m(k p d) S m[^-(x - w/2)} = M x (6.41) 

p w 

To find B p , we multiply both sides by sin [^(z - to/2)] and integrate over the extent 
of the aperture. As before, by employing the orthogonality relation 


t£ fl _ r , 

[ W/2 sin [— (z - to/2)] sin [— (z - w/2)]dx = j 2 
J-w/2 w w I o pi i q 


(6.42) 


we find 


2 n 


Q _ . 

p jki,w s'm(k p d) J-w/2 


f M x (x) sin [— (z — w/2)]dx (6.43) 

J-w/2 W 


The magnetic field in the internal region may now be expressed explicitly as 

*<*•»> = ^|^Mj c ° slMy + ' )|sinI - (l ‘“' /2)l 


f M r (z') sin [— (z' — w/2)]dx' 
J—w/2 tU 


which upon setting y = 0 gives the tangential aperture field 


H x (x, 0) = 


sin [— (z — to/2)] 


Yb 00 

kb tan (k p d) 1 to 


(6.44) 


r* MJx') sin [^(x' - w/2)]dx' 
J-w/2 tO 


Comparing this with (6.34) we deduce that 


(6.45) 


G e (z;z') = -E 


2k„ 


sin [— (z — to/2)] sin [— (z' — to/2)] (6.46) 


w tam(k p d) 1 w 


r p* 

u; 
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Substituting (6.32)-(6.34) and (6.46) into (6.3) yields the integral equation 
sin = Ml + - x'|)<ix' 

+ 1 f"' 2 sin p(x' - tx/2)]dx' (6.47) 

WTJ ej> J-w/2 W 


where 


rjep = JT~ z b tan (k p d) (6.48) 

K p 

are the normalized E-mode impedances of the cavity. Equation (6.47) is an exact 
integral equation to be solved for M x {x). 

Clearly, (6.30) and (6.47) are both invalid when 


tan(Ayf) = 0 (6.49) 

and this occurs only when the material filling the groove is lossless. To be specific, 
the modal solutions fail if there exist integers p and q such that 1 

( £ )! + ( j )i = ( f )J p ’ ,e:r (6 - 50> 

This difficulty in the evaluation of the internal Green’s functions may be circum- 
vented by assuming a small loss in the material. We also note that for the proper 
behavior of the field in the internal region, we must have 


&e{* p }> 0 (6.51) 

9m{i p } <0 (6.52) 

when usk-r (6.20). 

‘The formulation for the H-polarization case also fails if k p = 0 in addition to (6.49). This is 
equivalent to p/w = q/d when p,q £ 2. 
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Upon a solution of the integral equations, the scattering echo widths of the groove 
may be cal culated from (3.44) and (3.45) with Zjifj and Kx replaced by and 

2 M x , respectively, and the polarization subscripts e and h interchanged in accordance 
with the image theory and Babinet’s principle [19]. Thus, 


<?h = K 


/ w/2 
■w/2 


e }k a x' co« <f^ x ' 


(6.53) 


a. = k 0 


sin 


<t> r' 2 M x {x)e? koXf dx 

J —w/2 


(6.54) 


6.2.3 Numerical Solution 


The integral equations derived in the previous subsections may be solved numer- 
ically by the moment method and will serve as the reference for the validation of the 
results obtained from alternative formulations presented in the rest of this chapter. 

Considering the H-polarization case, the integral equation (6.30) may be dis- 
cretized by expanding M z (x) as 

M z (x) = ^2 M z (x n )P(x - x n ), x n = nA + — (6.55) 

n=0 

where P(x) denotes piecewise constant basis function. Substituting for the current 
expansion in (6.30) and applying point matching, the admittance elements are given 
by 


Imn — T mn II r 


(6.56) 


where T mn are elements attributed to the external tangential fields and are given by 
(3.85), while II mn are those attributed to the internal tangential fields and are given 
by 


n m n - A. £ — 1— cos p( t . - W2)] COS [£(*. - u,/2)]sinc(£^) (6.57) 
£ p wk p r) hp w w w l 
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A similar discretization can be carried out for E-polarization. 

Figure 6.3 shows a sample calculation of the backscattering echo width for an 
empty groove based on the above formulation. The groove is assumed to be 10A 
long and the usual physical optics approximation (1.26) was invoked to relate the 
two-dimensional echo width given by (6.53) or (6.54) to the corresponding three- 
dimensional radar cross section. The results are in good agreement with a corre- 
sponding finite element method (FEM) solution [67]. 

6.3 Partially Loaded Grooves 

We now consider the partially loaded groove shown in Figure 6.4. If the filling 
material is electrically dense, we may consider the equivalent problem of a homoge- 
neously filled groove of depth d terminated with a floor consisting of an impedance 
sheet. In this case, the boundary conditions on the cavity walls and the aperture 
of the groove remain the same as (6.10) and (6.11), while the floor satisfies the 
impedance boundary condition (5.58) 

E - (y • E)y = r\ t Z 0 y x H (6.58) 

where t)t is the normalized surface impedance of the floor. The above condition 
replaces (6.9), and in scalar form 

Ex = T} t Z 0 H z y — — d (6.59) 

E z = -veZoHr y = —d (6.60) 


Following an analysis similar to that of the previous section, new integral equa- 
tions can be derived for the solution . . ae equivalent magnetic current density over 

the aperture of the groove. In particular, employing the equivalence principle, the 
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Infinite 

ground 




Figure 6.3: Comparison of the backscattering patterns of a long two-dimensional 
groove obtained from a finite element solution (FEM) [67] and the 
method of moments (this study). The groove is assumed to be 10A long. 
(20 samples/ A with 60 waveguide modes). 
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Impedance Sheet, r| 
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Figure 6.4: Scattering from a groove partially loaded with electrically dense material, 
(a) Geometry, (b) Equivalent problem using an impedance sheet. 
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tangential magnetic fields in the external and internal regions are expressed in terms 
of the pertinent Green’s functions and used to enforce the continuity of the fields 
across the aperture. In this case, the expressions for the external fields remain un- 
changed and are given by (6.13) and (6.33) for H- and E- polarizations, respectively. 
As for the internal fields, the Green’s functions G h and G e must be modified to 
accommodate the new boundary condition on the groove s floor. 

Considering the H-polarization case, once again the wave potential ty TE is ex- 
pressed in terms of an infinite sum of orthogonal modes as in (6.16). The new 
boundary conditions to be satisfied by these modes are 

= r}tZ 0 (jhY b xl>$) y = -d (6.61) 

—tjj h =0 x = ±u>/2 (6.62) 

dx p 

A set of eigenfunctions satisfying the second condition along with the wave equation 
is 


= [e*M»+0 _ R h e~ ikpiv+d) ] cos [^j-(x - w/2)\ (6.63) 

where Rh is the reflection coefficient of the floor and k p is defined in (6.20). Enforcing 
the second condition yields 

jk,( 1 + R h ) = r,,J^-jk b (l - R„) (6 64) 

V Mr 


which upon solving for Rh gives 


Rh = 


Vi + ~r~ z f> 

k b 


(6.65) 


Upon imposing the equivalent current condition (6.11) on the aperture and solving 
for the new set of coefficients A p , we find 

2 j r w / 2 


A p — 


— rr f M z ( l)cOS [— (x — w/2)]dl 

£pwkp(e }k P d + Rhe~ }k i> d ) J-w/2 1 to V 


( 6 . 66 ) 
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and, therefore, the new Green’s function is given by 

G ‘ (x ’ 10 = § ^ [l£> + 008 1 - “Wl cos if-' 1 ' - u ’/ 2 )l 


(6.67) 


where e p is given by (6.29). Employing the above results, the integral equation to be 
satisfied by the aperture current for the H-polarization case is 


B jk 0 x cot <t > 0 


= ^ I”' 2 M,(x')H<?\k 0 \x - x'\)di 
4 J-wl 2 


cos [^0* - W2)] f w/2 


L XL± 

p=o £?Wh P 


J ^ M x (x') cos [~( x> ~ w/2))dx' (6.68) 


where rj' hp is the normalized equivalent surface impedance of the groove looking into 


the aperture, given by 


k v , , 

Vt + JZbT~ tan Kpd 

t Kh 

Vh P = * ^ * 

jr)t— tan k p d + z b 
k p 


(6.69) 


For the E-polarization case, the boundary conditions to be enforced in the internal 
region are 


V’p = ritZ 0 (jk b Y b -^tl>') y = -d 


v>: = 


x = ±w/2 


(6.70) 

(6.71) 


suggesting the following form for the eigenfunctions xfri 


tl>; = [c ifc » ,(v+<<) + R e e- jk *to + V] sin [^(x - u>/2)] 


(6.72) 


where R e is the reflection coefficient of the floor given by 


R c = 


k 6 

t)e - —z b 

~ *6 
( + T-Zb 

k„ 


(6.73) 
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Solving for the mode amplitudes B p , we find 

-2 Y b P/J 


B v = - 


- ~ ' rr- r r M x (x)sm[^-(x-w/2)]dx (6.74) 

jwkb(e ]kpd + Re e ,kpd ) J-w/2 w 


and the new Green’s function is given by 


<?•(*.*') = E 


-2k, 


P=1 




e^ kpd - R e e~ jkpd 


ei kpd + R e e * kpd 


sin [— (x — iu/2)] sin [— (x' — w/ 2)] 
w w 


rP 71 " / 


(6.75) 


Employing the above results, the integral equation for the E-polarization case is 
obtained as 

= t ( i + y&) - A)ix ‘ 

^ sin V^-lx - w/2)] H 2 / A . , /oxlJ , ^ 

4- — / M x (z)sin [ — (x — w/2)\dx (6.76) 

J—w/2 W 

where rj f hp is the normalized surface impedance of the groove looking into the aperture, 
given by 


rjt + 3 z bT- tan k p d 

f _ 

Vcp “ Z *> fc 

jytir tan ^d + z h 

kb 


(6.77) 


It is noted that (6.68) and (6.76) are identical to (6.30) and (6.47) with the only 

modification that the normalized mode impedances rjh p and T] ep are replaced by the 
normalized equivalent impedances Tj f hp and 7/' p looking into the aperture. It is also 
noted that if the groove is terminated by a perfect conductor, r\t — 0 and 


Vhp — Vhp j y ep — y*v 


and the formulation reduces to that of the homogeneously filled rectangular groove. 
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6.4 Dominant-Mode (Quasi- Static) Formulation 

For large apertures, a numerical approach is the only alternative to the solu- 
tion of (6.30) and (6.47). However, in many cases the characterization of narrow 
width grooves is of practical interest. With this motivation the narrow groove has 
been modelled as an impedance insert in an effort to simplify the analysis [68]. Un- 
fortunately, the resulting quasi-static integral equations were not amenable to an 
analytical solution but, nevertheless, it was possible to derive accurate empirical 
echo width formulae through the examination of numerical data. This was essen- 
tially done without a direct (analytical) evaluation of the current on the impedance 
insert. 

In this section we consider the solution of the integral equations for a narrow 
rectangular groove without invoking the impedance approximation used in [68]. It 
is shown that by retaining the dominant mode supported by the rectangular groove, 
the resulting quasi-static integral equations are comparable to those associated with 
the perfectly conducting narrow strip considered in Section ( 3.3.2). They are there- 
fore amenable to analytic solution yielding the exact field distribution or equivalent 
currents across the groove’s aperture. The derived currents exhibit the an edge be- 
havior similar to that associated with the currents of a perfectly conducting half 
plane or strip. On the other hand, the corresponding current behavior based on 
the (numerical) impedance simulation of the groove is quite different. However the 
resulting echo widths are comparable. 

The derived analytical expressions for the equivalent aperture currents are of 
potential importance in constructing suitable models for long and narrow three- 
dimensional apertures. A so, unlike the echo width formulae given in [68], th. se 
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derived here are valid for all groove depths and material fillings. In this sections, 
the exact integral equations derived for a two-dimensional rectangular groove will 
be simplified to the case of a narrow width groove and solved for the equivalent 
magnetic currents across the aperture. The accuracy of the currents is examined by 
a comparison with the numerical data. Simple echo width expressions are also given 
for the principal polarizations which are treated separately. 

6.4.1 H-polarization 

When kw Cl, the Green’s function G h (x; x') can be substantially simplified by 
retaining the first term of the sum, corresponding to the lowest order mode in the 
cavity. The integral equation (6.30) then reduces to 

= -J- r ,! r M,{x’)HV>(k.\x-x'\)ix' 

2Wf)h. J-w/2 4 J-w/2 

(6.78) 

where 

Vh = jJ^ta.n(k b d) (6.79) 

is the normalized impedance of the dominant mode which can be also identified as 
the normalized impedance at the surface of a grounded slab of thickness d. It is 
interesting to note that if the first integral in the right hand side of (6.78) is replaced 
by wM z (x ), then (6.78) reduces to 

z _ MM + h f" n - x'\ )dx' (6.80) 

2 T)h 4 J-w/2 

which is the integral equation based on the impedance boundary condition 


E x (x) = tj h Z 0 H z (x ) 


(6.81) 
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applied over the extent of the aperture (see equation (5.32)). On the other hand, the 
integral equation (6.78) is based on the relation 

1 rW 2 

- J E x (x')dx = r] h Z 0 H z (x) (6.82) 

This observation reveals the inherent local nature of the impedance boundary con- 
dition and its underlying assumption that the current is more likely to be slowly 
varying. Not surprisingly, (6.80) predicts a rather smooth behavior of the current 
distribution near the edges of the groove at x = ±w/2. In contrast, a numerical 
solution of (6.78) gives the usual singular form of M z (x) at the same locations. An- 
other interesting property of the “boundary condition” (6.82) is the independence of 
the left hand member from x. This property will be exploited later to arrive at a 
closed-form solution to the integral equation. 


To further simplify (6.78) for k a w C 1, we introduce the small argument expan- 


sion for the Hankel function as we did for the narrow strip, 


Hi 2 \z) ~ 1 “if In (2p) + C?(z 2 ,z 2 \nz) 


(6.83) 


where In 7 is Euler’s constant. Substituting this into (6.78) and retaining only terms 
to O(fcw) we have 

/ W 2 9 * 

M z (x f ) In |x — x'\dx' = j^—Z 0 
•w /2 k 0 

-'(¥)-<§]£>■•»■• <“> 

Further, by introducing the same change of variables as ( 3.51) 


_ 2x , _ 2x' 


(6.85) 


(6.84) becomes 


/-. b + fe - ,n (^f 2 ) - ^ /-’ 


( 6 . 86 ) 
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The above singu lar integral equation can be inverted by noticing that the right hand 
side of (6.86) is independent of £, and upon invoking the identity ( 3.54) 

J (1 — x' 2 ) _1 ^ 2 ln \x — x'\dx' = — 7rln2 x€[— 1,1] 


we 


find 


M,(f) = Wi-{T 1/2 = ^ 




(6.87) 


where Xh is a complex constant given by 


Xh = 


4; 




( 6 . 88 ) 


It is noted that the aperture magnetic current (6.87) has a functional form exactly 
similar to that of the electric current of a narrow strip ( 3.56). In fact, when rj h -* oo 
corresponding to an open slot, we find that 


lim Xh = 

fJh—oo 


4j 


(6.89) 


^R¥m 

which is analogous to the E-polarization result obtained for the narrow strip (3.58). 
This result is, of course, expected based on Babinet’s principle [19, 66]. 


6.4.2 El-polarization 

A si mil ar derivation can be carried out for E-polarization. By retaining the lowest 
order mode in the groove, (6.47) becomes 

sin4> 0 e’ k * am *° = — — cos(— ) f M x {x')cos dx' 

WT] e W J-w/2 W 


4 


(6.90) 
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where 

tan(Arid) (6.91) 

is the normalized wave impedance of the lowest order mode with 

*?- **"(£)' (6 - 92) 

Next, by introducing the small argument expansion of the Hankel function and the 
change of variables defined in (6.85) we obtain 

sin A, = ± C0S (| f )/-, M-M 

(693) 

where we have retained only terms to O(kw). An approximate solution for M x can 
now be obtained by satisfying (6.93) at £ = 0. We have 

^ Li Mx ^ ln K “ = j*k 0 w sin j_ x Af x (£') cos(|^')^' 
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and J\ denotes the first order Bessel function. In arriving at this result, we also 
made use of the identity [69] 

/ \/l — x 2 cos (nx)dx = —Ji(n) (6.97) 

J - 1 W 

It is noted that r) e is generally inductive for kw < 1, and thus the derived 
expression for Xe h nonsingular within the expected validity range of (6.95). Again, 
as T) e — ► oo corresponding to an open gap, 


lim \h = jk D w sin <f> 0 (6.98) 

TJc-^OO 


which is analogous to the H-polarization result obtained for the narrow strip ( 3.59). 

The far zone scattered field at a point (p, <j>) in cylindrical coordinates can be 
computed from (6.53)-(6.54). Upon approximating the exponential e■' fcoX ' co, * with 
unity, we have 




2 

, 7TW 

<r* = k 0 

Y 0 / M x (x )dx 

J-w/2 

= k 0 \—Xh 


(6.99) 


and 


I rwfl 

sin <f> / M x (x')dx' 
J—xv/2 

which in the backscattering direction, yield 


= k. 


icw 


Xe sin 4>\ 


( 6 . 100 ) 


2xA 


CTh = 


*'§(■*&) 


( 6 . 101 ) 


7T A 


= 


( k 0 w y 

\2sin^/ 


2 1 + 0.5696 


jk 0 w 


( 6 . 102 ) 


Ve 


Before a detailed examination of the the above quasi-static results, we remark 


that the same analysis presented above is applicable to a narrow groove whose floor 


Figure 6.5: Geometries of some gaps and crack of practical interest. 

satisfies an impedance boundary condition. In this case the mode impedances and 
Tj e are replaced by the corresponding normalized equivalent impedances Tj‘ h and r/' 
looking into the aperture. This allows an analysis of partially filled narrow grooves 
as well as narrow cracks of simple shapes (Figure 6.5). For such geometries, a quasi- 
static or empirical estimate of the impedance may be used [68]. 

The derived formulae for the gap echo width and aperture currents are based on 
low frequency approximations to the exact integrals. They are thus expected to be 
valid for small groove widths and it is, therefore, of interest to examine their accuracy 
limitations as the width of the groove increases. Also, of interest is a comparison 
of the analytical echo width formulae derived here with the corresponding empirical 
ones given in [68]. 

Figure 6.6 presents a comparison of the derived H-polarization current distribu- 
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tion (6.87) versus that obtained from a numerical solution of the full-wave integral 
equation (6.80). Similar comparisons are also given in Figure 6.7 for E-polarization. 
In both cases <f> 0 = ir/2 and for this incidence the expressions (6.87) and (6.95) are 
in good agreement with the exact data (although only amplitude comparisons are 
shown, good agreement was observed for the phase as well). This holds independent 
of <t> 0 for small w. As the groove width increases, however, the exact current is to 
an increasing extent a function of <j> 0 and as noted in [68] the angular dependence 
is noticeable for w > 0.15 A. Since the quasi-static H-polarization current (6.87) 
is independent of <t > 0 , it is then applicable up to this value of w. Nevertheless, we 
have found that for normal incidence, (6.87) is quite accurate up to w « 0.25A and 
its accuracy improves for filled grooves. For E-polarization, the derived quasi-static 
current solution is an explicit function of <j> 0 and, therefore, remains accurate for all 
angles of incidence up to w 0.25A. 

Comparisons of the echo width formulae with numerical data are given in Figures 
6.8 — 6.11 for the H-polarization and Figures 6.12 — 6.13 for the E- polarization case, 
respectively. These results correspond to the backscattering computations at normal 
incidence (<^> = <f> 0 = 7t/2). It is observed that the quasi-static formulae remain 
accurate for all groove depths provided w is kept within its validity bounds. The 
empirical formulae given in [68] were generally found to agree with these results, 
except near the resonance regions for the H-polarization where the empirical formula 
fails. This is illustrated in Figure 6.11 for an empty groove whose resonant depth is 
d = 0.234A when w = 0.1A. Also, in contrast to H-polarization, the E- polarization 
echo width does not display any resonant characteristics for small w since there is 
no traveling mode in the cavity. In fact, for w < 0.2A, the E-polarization echo width 
of an empty groove is independent of depth for d > 0.1 A. 
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Finally, we remark that the above solutions are of potential utility in the analysis 
of long three-dimensional (finite) grooves. For example, Figure 6.14 shows the radar 
cross section from a 2.5A long groove whose width and depth are A/4. In this case, 
the quasi-static result was obtained from (1.26) based on the physical optics ap- 
proximation. Good agreement with the full-wave three-dimensional moment method 
solution (Chapter 7) is observed. 
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Figure 6.6: H-polarization equivalent surface magnetic currents for a groove of width 
w = 0.1 A and depth d = 0.2A; Comparison of analytical and numerical 
data. 



Magnetic Current |Mx|. A 


164 
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Figure 6.7: E-polarization equivalent surface magnetic currents for a groove of width 
w = 0.1 A and depth d = 0.2A; Comparison of analytical and numerical 
data. 
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Figure 6.8: H-polarization normal incidence echo width for a groove of depth d — 
0.2 A as a function of width for three different material fillings. 
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Figure 6-9: H-polarization normal incidence echo width for an empty groove as a 
function of depth for three different widths ( tv = 0.05A,0.1A, and 0.2A). 
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Figure 6.10: H-polarization normal incidence echo width as a function of depth for 
a groove of width w = 0.15A filled with a material having t r = 4 — jl 
and n T = 1.5 — j’0.05. 
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Figure 6.11: Comparison of the quasi-static and empirical solutions [68] with numer- 
ical data as a function of width for d = 0.234A (near resonance). 
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Figure 6.12: E-polarization normal incidence echo widths for an empty groove as a 
function of width (<f = 0.2A); Comparison of the quasi-static, empirical 
[68], and numerical solutions. 
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Figure 6.13: E- polarization normal incidence echo widths for an empty groove as a 
function of depth; Comparison of the quasi-static, empirical [68], and 
numerical solutions. 
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Figure 6.14: Backscattering cross section of a long empty cavity (£ = 2.5A, u; = 
0.25A, d = 0.25A); Comparison of the quasi-static (factored echo width) 
and full-wave three-dimensional Moment Method solution. 
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6.5 GIBC Formulation 

In the previous sections, we presented a rigorous full-wave formulation for com- 
puting the scattering by a filled rectangular groove in a ground plane. This was 
further approximated to the case of a narrow groove based on a quasi-static anal- 
ysis of the pertinent integral equations. In this section, we present another class 
of approximate formulations for the general analysis of grooves which make use of 
generalized impedance boundary conditions. 

As mentioned in Chapter 5, the material filled groove may be simulated by a 
two-dimensional impedance insert. Indeed, we have already encountered the SIBC 
formulation for the groove in equation (6.80). In this section, we examine the accu- 
racy of this boundary condition as well as those of higher orders. 

Consider first the SIBC. In this case, the integral equations (5.32) and (5.33) axe 
applicable and by setting F t = F' t = 0, l = 2,3 in (5.34) and (5.35), we obtain 

Z 0 t ik ***~+* = ^Fi(x)M,(x) + DFT -1 {m,T} (6.103) 

and 

•m^ 4 — *• = + jj DFT "> {(*’ - 

(6.104) 

which are subsequently solved by the CGFFT method. 

Figure 6.15 shows the amplitude and phase of the equivalent magnetic current 
density for a half-wavelength deep (in free space), two wavelengths wide rectangular 
groove filled with a lossy material of high index of refraction. In this case, the 
conditions (5.5)-(5.6) are satisfied 


|N| = 7.7 |3mN|M = 5.5 
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and the agreement with the full- wave solution is very good. In contrast, when the 
groove is filled with a relatively low contrast material, SIBC is no longer applicable 
as the validity conditions of SIBC are violated (Figure 6.16). 

Consider the same groove, now simulated by an impedance sheet (insert) satis- 
fying a third order boundary condition. In this case, the integral equations (5.23) 
and (5.31) are applicable and the CGFFT implementation is given by (5.34)-(5.35). 
From Figure 6.17, it is seen that the GIBC solution agrees reasonably well with the 
exact one except near the groove terminations. 

Generally, the current distribution based on the proposed third-order GIBC is 
not of acceptable accuracy when within 0.25A of the groove’s terminations. However, 
because it is in good agreement with the exact current distribution elsewhere, one 
approach in retaining the memory advantage associated with the application of the 
GIBCs is to combine the exact and GIBC formulations. This is discussed in the next 
section. 

6.6 Hybrid Exact-GIBC Formulation 

Based on the above discussion, a procedure for combining the exact and GIBC 
formulations is to feed the currents predicted by the GIBC integral equation (5.23)- 
(5.31) away from the edges into the exact integral equation (6.30)-(6.47). The last 
can then be solved for the remaining currents in the vicinity of the groove termi- 
nations. This only requires the inversion of a small matrix and hence the memory 
demand is essentially unaffected. 

To demonstrate this hybrid approach, let us consider the H-polarization and a 
similar formulation applies to the other polarization as well. Suppose that Mf(x) 
denotes the current computed via the GIBC integral equation (5.23) and we choose 
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Figure 6.1 5i Simulation of a groove filled with a high-contrast material; w = 2 A, 
d = 0.5A; e r = 12.5 — j2.5, fi r = 4.5 — jl.2, and <f> 0 = 30°. Comparison 
of the full- wave ( — ) and SIBC ( ) results. 






A , dB 













Figure 6.17: Comparison of the full- wave (— ), SIBC (- - -), and the third order GIBC 
(• • •) solutions for a filled groove; w = 2A, d = 0.5A; e r = 2.5 - >0.2, 

fir = 1.8 - >0.1. 
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to approximate the true aperture current as 


fMf(x) |x| < f - x A 


M,{x) = { 


(6.105) 


{ Mf (x) |x|>*-x A 

where denotes the unknown currents near the edges of the groove. To compute 
we substitute (6.105) into (6.30) and this yields 


2e , *'* c "*"+ M°( x') [;'i„HGi(x,x') - ^H^>(k 0 \x - x'|)j dx' 

= ~ Mf(x') [jhY t a h (x,x‘) - ^H^(K\x - x'|)] dx' 

- A f,V) [i*»HG»(x, x') - *f*Hj J) (*;„|x - x'|)] ix' 

(6.106) 


Assuming that M z G (x) has already been computed via a CGFFT solution of (5.23), 
the entire left-hand side of (6.106) is known and thus, for xa < 0.25, a 4 x 4 or a 
6x6 square admittance matrix is required for the solution of M G (x). In general, 
continuity of the current density must be imposed at the transition regions between 
Mf(x) and Af G (x), and this can be accomplished through a simple averaging. 

Figure 6.18 shows the results obtained for the aperture current density of the 
groove considered before. Clearly, the proposed hybrid solution (HYBRID-3) pro- 
vides the necessary correction near the terminations where the GIBC solution fails. 
Bistatic and backscattering patterns corresponding to the same groove are given in 
Figure 6.19. It is observed from these patterns that the SIBC solution is substan- 
tially in error for angles near grazing. The same holds for the GIBC since, as is well 
known, the contribution of the edge currents is a dominating factor in the compu- 
tation of the echo width. Notably, the patterns predicted by the hybrid formulation 
are always in good agreement with the full-wave moment method solution. 
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Figure 6.18: Comparison of the current density for a rectangular groove as obtained 
by the full-wave and approximate formulations; w = 2A, d = 0.5A; 
^ = 2.5 — jO.2, fi r = 1.8 — jO.l, and 4> 0 = 30°. (a) Amplitude, (b) 
Phase. 


c 
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Figure 6.19: Comparison of TE scattering patterns for a rectangular groove as 
obtained by the full- wave and approximate formulations; w = 2A, 
d = 0.5A; t T = 2.5 - >0.2, fi r = 1.8 - j'0.1. (a) Bistatic <f> 0 = 30°. 
(b) Backscatter. 
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6.7 Tapered Grooves 

The GIBC formulation can be directly applied to the scattering from tapered 
grooves, provided the constant coefficients associated by the employed GIBC are 
allowed to vary. This, clearly, avoids a need to compute the Green’s function or to 
use a more sophisticated technique such as the finite element method (FEM) [70]. 
The condition on the slow variation of the impedance for the SIBC is given by (5.7) in 
addition to (5.5) and (5.6). However, it is possible to simulate more rapid variations 
by using & higher order GIBC. Consider, for example, the non-rect angular groove 
shown in Figure 6.20. In this case, the SIBC is inadequate in modelling the groove 
while a direct application of the third order GIBC formulation is sufficient to yield 
accurate results. 

6.8 Summary 

The problem of scattering from two-dimensional rectangular grooves was studied 
using a full-wave analysis. The analysis is applicable to grooves terminated with 
perfect or imperfect surfaces. This formulation was specialized to the case of elec- 
trically narrow grooves by considering the dominant waveguide modes in the groove 
and employing the finite Hilbert transform theory based on a quasi-static approxima- 
tion of the resulting integral equations. Analytical expressions were derived for the 
equivalent magnetic current distribution over the aperture of narrow grooves. The 
solutions were found to exhibit the familiar edge behavior observed in the case of 
narrow strips and slits. Using the derived current distributions, closed form expres- 
sions were given for the echo width of the narrow groove and these were compared 
with numerical data. Their accuracy was examined as a function of width, depth 
and material filling and were found to be in good agreement with the echo width 
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w-lA. WJ.03A, VA-055, H-Pol. 



Angle of Incidence deg. 



Angle of Incidence p , deg. 


Figure 6.20: Comparison of TE backscattering results based on the SIBC, GIBC-3, 
and FEM [70] for a non- rectangular filled groove (e r = 5 ,/x r = 1). (a) 
Geometry, (b) w = lA,d = 0.03A,6 = 0.25A. (c) w = 2A ,d = 0.1 A, b = 
0.5A. 
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data based on full- wave solution for all angles of incidence, provided w < 0.15 A for 
H-polarization and w < 0.25A for E-polarization regardless of the groove’s depth. 
The closed form solutions obtained here were found to be of potential use in the 
study of the long and narrow grooves and could significantly simplify their analyses. 

Furthermore, the scattering behavior of the groove was simulated by the impedance 
boundary conditions. Both first order (SIBC) and third order GIBC formulations 
were studied. The formulations based on these boundary conditions were found easier 
to implement than the full-wave formulation. Also, unlike the exact integral equa- 
tions, they were amenable to a CGFFT implementation. For high-contrast material 
fillings, the SIBC was found adequate in modeling the groove. An analytical compar- 
ison of the integral equation based on a SIBC simulation with that from a full-wave 
formulation, revealed a well-known limitation of the SIBC formulation. That is, the 
SIBC integral equation generates an average of the actual current distribution. By 
resorting, though, to a third order GIBC the correct current behavior was reasonably 
predicted away from any abrupt terminations of the groove. The predicted current 
based on the GIBC simulation was in general incorrect near the edges and to correct 
this deficiency, a hybrid approach was proposed. Specifically, the currents computed 
via the GIBC formulation away from the rectangular groove terminations were em- 
ployed in the exact integral equation to generate a small matrix for the currents 
in the vicinity of the terminations. This was referred to as the hybrid exact-GIBC 
formulation and was found to yield a reasonably good prediction of the scattering 
by filled rectangular grooves. 

Finally, when the groove terminations are not abrupt, the hybrid formulation is 
not required and a direct application of the GIBC formulation may be sufficient. 



CHAPTER VII 


SCATTERING BY OPEN RECTANGULAR 
CAVITIES RECESSED IN GROUND PLANES 


7.1 Introduction 

The characterization of apertures in a ground plane is of considerable importance 
in radar cross section (RCS) and electromagnetic pulse (EMP) studies. Indeed, a 
large body of work exists for the analysis of two-dimensional slits in a thick ground 
plane [71]-[76] or cavity-backed apertures [77]-[79]. Extensions of these procedures to 
three dimensional characterizations are possible, but so far this has only been done 
for high frequency techniques. Numerical solutions for three dimensional apertures 
have been limited to scattering and transmission by openings in a thin ground plane 
[50], [80]- [82] primarily due to the excessive computational demands and complexity 
of the solution. The only exception to this is the use of the mode-matching tech- 
nique for the analysis of rectangular [67] and spherical [83] cavity-backed apertures. 
Although in principle exact, the mode-matching approach leads to an infinite system 
of equations in addition to being cumbersome. A need, therefore, exists to develop 
numerical solutions for cavity backed apertures. Such solutions can provide a charac- 
terization of this structure and could serve as a reference for validating new solution 
algorithms. 
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In this chapter we consider the scattering by a rectangular cavity-backed aperture. 
The solution technique employed in the analysis is the full-wave moment method ap- 
proach considered in the two-dimensional applications of Chapter 6. A fundamental 
aspect of this method is to employ the aperture fields as the equivalent sources of 
the fields interior and exterior to the cavity. The complete integral representation of 
the fields within the cavity makes use of the modal Green’s function whereas that 
external to the cavity makes use of the free space dyadic Green’s function. An inte- 
gral equation for the aperture fields is then constructed by enforcing tangential field 
continuity across the aperture. Except for being tedious, the entire solution process 
is straightforward and in an effort to maintain the level of complexity to a minimum, 
a pulse-basis moment method solution of the integral equation is first discussed. The 
more useful roof-top basis is also presented. As can be expected, the roof-top basis 
formulation leads to a more efficient numerical solution at the expense of additional 
complexity. In either case, the admittance elements associated with the external 
fields are identical to the impedance elements for a perfectly conducting plate. How- 
ever, the major difference in computational efficiency among the two formulations 
lies in the evaluation of the admittance elements associated with the internal fields. 
These are given in terms of a double sum series whose convergence is substantially 
improved when higher order basis functions are employed. 

In the following sections we first develop the complete field representations in the 
interior and exterior regions of the cavity. The integral equation is then formulated 
by requiring continuity of the tangential magnetic field across the aperture and dis- 
cretized using pulse and roof-top basis functions. The evaluation of the admittance 
elements is discussed in some detail since these are of crucial importance in the over- 
all accuracy and efficiency of the solution. An important aspect of this chapter is 
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Figure 7.1: Geometry of an open cavity recessed in a ground plane. 


the presentation of a number of scattering patterns some of which are validated with 
data obtained via an alternative solution method. 

7.2 Full- Wave Formulation 

Consider the aperture shown in Figure 7.1 illuminated by a harmonic plane wave 
given by (5.36) and (5.37). This represents a two-media problem with the aperture 
dividing the space into two regions, one external to the cavity (z > 0) and another in- 
ternal to it (— c < z < 0). To formulate the fields scattered by the cavity, an analysis 
similar to the two-dimensional case is carried out based on the equivalence princi- 
ple. Accordingly, the aperture is closed by a perfect conductor and the equivalent 
magnetic current 

M = Exn = Exz = xE y — yE x (7.1) 

is placed on the aperture at z = 0 + . The radiation of this current represents the 
scattered field in the external region and by demanding continuity of the tangential 
electric field, it follows that equivalent sources for the internal fields are the magnetic 
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currents 


M' = E x n' = E x -z = -M (7.2) 

placed across and just below the aperture at z = 0 - . It remains to also enforce 
continuity of the tangential magnetic field across the aperture and this will provide 
the required condition for determining the magnetic currents. By denoting the fields 
in the external region as (E°,H°) and those in the internal region as (E\ H 6 ), we 
have 


z x (H°(M) + H‘°] = z x H fc (M'), z = 0 (7.3) 


where H ,a represents the incident magnetic field in the absence of the aperture. Upon 
substitution for H° and H\ (7.3) then yields an integral equation for the equivalent 
magnetic currents. 

The external scattered field can be expressed as the surface integral 

H°(r) = -jk 0 Y 0 [ 2M(r') • f (r; r')ds' (7.4) 

J s 

where S denotes the surface of the aperture, T is the free space dyadic Green’s 
function 


f(r;r')= (l + ^wW(r;r') 

Also, a factor of two was introduced in (7.4) to account for the ground plane’s 
presence. By expanding the T in cartesian coordinates, (7.4) can be written more 
explicitly as given by 


HI 

HI 


2jY 0 

K 

2jY 0 

k 0 




G 0 (r;r')ds' (7.5) 
G 0 (r;r')ds' (7.6) 
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The internal fields can be written in terms of the TMzand TEz waveguide modes. 
We have 


E 6 = E™ + E TE = -jk b Z b 


1 +XF VV. 


(z^™) — V x (£^ rE ) 


1 


1 + I vv. 


(z'V**') 


(7.7) 


(7.S) 


H 6 = H™ + H te = V x (z9™) - jk b Y b 

where as before, k h = Kk 0 is the wave number in the internal region and Z b * 1/ft 
is the intrinsic impedance of that region. The functions ¥™ and <H TE are the wave 
potentials both satisfying the wave equation 

(£*£*&*■>)— 

subject to the boundary conditions 


(7.9) 


E x = E y = 0, 2 = -c 

(7.10) 

E x = E z = 0, y = 0 and y = b 

(7.11) 

E v = E z = 0, x = 0 and x = a 

(7.12) 


on the cavity walls. Referring to Figure 7.1, we have 

V™ = f) f; A m „ sin (™x) sin cos [k z (z + c)} 

m=l n=l 

and 

= E E ft™ cos (^x) cos (y ») sin [Me + <=)] 


(7.13) 


m=0 n=0 


(7.14) 


excluding m = n = 0 

where A mn and B mn are constants to be determined and k z satisfies the separation 
parameter equation 

+ (^) 2 + k\ = ki (7.i5) 
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Observing the restriction on the mode propagation constant 7 = jk z for the proper 
field behavior, we demand that 


Re{fc z ) > 0 

3m{fc*} < 0 


(7.16) 

(7.17) 


implying that 


k g — k mn — \ 




mr\ 2 

T) 


. y ■ W(™)Mt ) 2 

when kb is real. Substituting (7.13)-(7.14) into (7.7) we obtain 

* - £ h (x) k ~f A - + (t) H 

cos ("£T X ) sin ("F y ) sin + 


and 


sin (“7" x ) cos (l T y ) sin ^ mn ( 2 + 


Hi = 


eKt 

sin cos cos [k mn (z + c)] 


008 ("a~ X ) Sin (T y ) C ° S ^ mn ^ 2 + C )1 


(7.18) 


(7.19) 


(7.20) 


(7.21) 
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The mode coefficients A mn and B mn can be expressed in terms of the equivalent 
current M by enforcing the boundary conditions (7.1) and (7.2) 


E x = M' = -M, 


v> 


2 — 0 


E v = -M' = M x , 2 = 0 


and by invoking mode orthogonality we find 


A mn — 


2jhYb ( ( 


mir \ 2 /titiA 2 ) 

—) + (t) } 


-1 


kmn ab sin (k mn c) ^ V 


Bmn — 




_i u 

ab sin (A? mn c) ^ \ 


-l 


where 


jmn 


jmn 


= jf M s (x', y') sin cos (“^V) ds> 

= J M y {x', y') cos (-jp-*') sin (-yV) ^ 


and 


£m = < 


1 m = 0 

2 m > 1 

The above expressions for the internal cavity fields are invalid when 

kmn tan (fcmnC) = 0 


(7.22) 

(7.23) 


(7.24) 


(7.25) 


(7.26) 

(7.27) 


(7.28) 


(7.29) 


which may occur if the cavity is filled with lossless material. Hence, the modal 
solution fails if there exist integers m, n, and p such that 




(7.30) 
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where N is the index of refraction of the material filling the cavity. As mentioned 
before in connection with the two-dimensional grooves, this situation may be handled 
simply by introducing a small loss in the material. 

The desired system of integral equations is now obtained based on the continu- 
ity of the aperture magnetic fields (7.3), upon substitution for the pertinent field 
quantities. 


7.2.1 Reduction to the Two-Dimensional Case 

Before we consider the numerical solution of the above integral equations, we 
consider the special case of the long cavity. We show that the above formulation 
reduces to that of the two-dimensional case studied in Chapter 5 when the cavity is 
taken to be infinite along one aperture dimension. Although this analysis is quite 
general in its application, we limit our attention to the H-polarization case here. In 
particular, as b -> oo, we may neglect the contribution of the transverse x-component 
of the equivalent magnetic current density in favor of the dominant longitudinal y- 
component. Also, the problem is invariant in the y-direction and we set — = 0. 

dy 

Hence, from (7.6) and (7.21) we obtain 

H > = -VWo J _ w/2 M y (x', y')G p (x, y\ x, y')dx'dy' (7.31 ) 

and 


lim Hi 

fc—oo v 


= -2 jk b Y b 


a b H Yh k 

uu m= 0 n= 0 "'mn 






jmn 


(7.32) 


Using the identity 


/“ + V r )dy = 


(7.33) 
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equation (7.31) can be rewritten as 

k 0 Y 0 r / 2 


lim HI = of, 

6—* oo ^ 2 J — 1 1//2 


f W/2 M y {x')H (2) {k 0 \x - x'\)dx' 
J — ur/2 


(7.34) 


which is compatible with the two-dimensional scattering integral (6.13). Further- 
more, substituting for J y mn and setting fc m „ = k m ^ Jkl - (*f ) 2 in (7.32) yields 

5 cos( HE x ) !' My{x ' )m< ™x')dx> (7.35) 
v a Jo a 


i; rri _ -j^kYb y' £m ,_ 

a “*0 fcm tan(fc m d) 


mil 11 J, 

6— *oo v 

where 


2 ~ . nir . y 00 . ,n7r , , 

S,-- b I>^T y) L 3m( T v)dv 

n=0 


(7.36) 


and S v should equal unity for (7.35) to reduce to the two-dimensional result (6.30). 
This can be verified by invoking the distribution theory. Specifically, by interchanging 
the order of summation and integration, and noting that [84] 


oo 

sin(mry) sin(n7ry') = 6(y — y') 

n=l 


(7.37) 


we find that 


= *>(?») ■fot'fy')] dy ' = jT <(» - = i (7.38) 


7.2.2 Numerical Solution via Galerkin’s Method 

In accordance with the method of moments, the integral equation to be solved 
for M is (see (7.3)) 

^[H?(x,y)-Hj(x,y)] •W(x,y)ds = -^H‘°(x,y)-W(x,y)ds , z = 0 

(7.39) 


where H? = x x H°, H b t = z x H\ and H|“ = z x H io with H° and H k as given by 
(7.5)-(7.6) and (7.20)-(7.21), respectively and W(x,y) is a weighting function. To 
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discretize (7.39), M is expanded in terms of subdomain basis functions 

M (*> y) = J2 m pi 4 ~ x p\y~ y,) + y - y,)] (7.40) 

where £ x and £„ are separable functions of x and y representing the expansion func- 
tions in the x and y directions, respectively. Further, M p , = xM xpq + yM vpq are the 
unknown coefficients of the basis functions. In accordance with Galerkin’s method 
we will choose the weighting functions to be the same as the expansion functions, 
i.e. 

r 

x( x (x - x,; y - y , ) 

W (*,y)=< (7.41) 

y( v { x x i,y i ij) 

for testing at the point (x,-,j/j). 
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where the admittance elements with the superscript a are associated with the external 
fields whereas those with the superscript b are associated with the internal fields. The 
external admittance elements Y* x , Y y a x and Y y a y are given in terms of the integral 

9H„ = s(*i, 

(7.45) 


= / ,ai [“* r ’ G a (r,r')dy'dx'dydi 

J{i-l)Ax J(j- l)Ay J(p-l)&x J(q-l)Ay 

and its second derivatives. This requires analytical evaluation when | i - p| < 1 and 

| j — g| < 1 and to do this, we rewrite gijpq as (i2 |r r |) 

r e -jk 0 R i 


9ijpq 


=-/ / 

4 7T J Sij w Sp<j 


ds'd, + -5- / [ ^ds'ds 

4tt Js„ Js„ ti 


R R\ 

The first integral has a nonsingular integrand for all i,j,p and q and can, therefore, be 
evaluated numerically using, for example, Gaussian integration. The second integral 
has a singular integrand when i = p and j = q but can be evaluated analytically to 


yield 



‘ (x -s'Xy-yQ 
2 


{(jr - x') In [( y - y ') + H] 


+(?/ - y) In [(x - x ) + R]} 


(x - x')(y — y')[(s - x') + (y - y')] 

4 


(7.46) 


^31 i pAx |?A» |«Ax tjAy 

_ g I r'=(p— l)Ax I y'=( 9 — l)Ay !*=(«— l)Ax I y=(j— l)Ay 

Unfortunately, the derivatives of gij pq do not exist in analytical form. A possible 

alternative is to evaluate them discretely using the computed values of gijpq, and a 

convenient way to do this is to employ the discrete Fourier transform. Proceeding in 

this manner, we first define the sample train 

N. N, 

xn = YdHsupqtii* - x p)Hy - Vi) 

p=l g=l 


( 7 . 47 ) 
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whose two-dimensional DFT will be denoted by Xij. Using a central difference scheme 
for the derivatives, the DFT of the sample train 

—~\ll 92 e ~ jkoR 

p=i 7=1 


N* r t r 

= l l 


- ds'ds 


->tj ^X 2 47r/? 

with x, and yj kept constant, can be approximated 


S(x - x p )S(y - y q ) 


(7-48) 


as 


Xij D'Xij (7.49) 

The admittance matrix elements can then be expressed as 


iXLh = - 

DFT -{[*>_ BJ] 

(7.50) 

fe), = 

2 iY 

3 0 DFT - 1 {-D x D y Xij} 

(7.51) 

II 

d 

2j Y ° DFT - {-D,D,x.,) 

(7.52) 

te), = - 


(7.53) 


where DFT -1 denotes the inverse discrete Fourier transform. These give the ad- 
mittance elements for the matrix row associated with the testing point The 

other row elements can be obtained by a simple rearrangement of this row upon 
invoking the symmetry properties of the matrix. 

To evaluate the admittance matrix elements associated with the cavity region we 

refer to (7.20)-(7.21). Substituting the expansion (7.42) into (7.24)-(7.25) and then 
into (7.39) yields 

(7.54) 

sm (— *,) sin cos cos 
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K) (J „ = -c EEw(=)(t) 

J ri m n 

/ TtlTT \ ( TTl 7T \ . / TUT \ f T17T \ 

■ cos ( v 1 ') sm v ~ Xi ) sm It * 1 ') cos ) 

• >in ( x 1 -) 108 ( v *') 1:08 ( t s ’) sin ( T y ’ ) 

and 

In these 

sine 2 Ax) sine 2 (%Ay) 

^ ^mn tan (fc TOn c) 

c _ 2jYo (AxAy) 2 
fc 0 /4 <*& 


( 7 . 55 ) 


( 7 . 56 ) 


( 7 . 57 ) 


( 7 . 58 ) 

( 7 . 59 ) 


and £ m have been defined in (7.28). 

It remains to compute the excitation elements and (/y) „ given by 

( v ttAx fjAy x . 

7“ c ) =2/ / H' x (x,y,z = 0)dydx 

rt Ax fjAy 

fr nc ) .. = 2 1 / #;(*, y, * = 0)<fydx. 

' V '*J ./(»— l)Ai J(j— l)Ay 

Integrating, we obtain 

_ 2/7 0I e jfco “ neo ^* iCO ** <,+Vj " n * o ^AxAy 


sine f 


k a Ax sin 6 0 cos <f>, 


sine ^ 


kg Ay sin fl 0 sin <ft 0 
2 


( 7 . 60 ) 

( 7 . 61 ) 


( 7 . 62 ) 
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and 



2H 0v e jko,iae °l Xi “•*•+» 


sine 


fc 5 Ai sin 0 O cos <j>, 


sine ^ 


(7.63) 


k 0 Ax sin 0 O sin <t> c 


where H ox and H oy are given in (5.38). 

This completes the derivation of all elements appearing in the system (7.44) 
whose solution yields the current densities M x and M y . Unfortunately, the compu- 
tation of the matrix elements Y xx , Y xy , Y yy and Y yx requires the evaluation of double 
infinite summations which are slowly converging. The asymptotic behaviors of the 

summation terms for Y xx and K 6 „ are of the form , = and , = 

y nVn 2 + m 2 mnv/n 2 + m 2 ’ 

respectively, implying that for fixed n, Y xy and Y yx have slow convergence whereas 
Y xx and Y yy are strictly non-convergent. As will be seen later, however, for long 
and narrow cavities, a finite number of summation terms are sufficient to obtain 
acceptable results. Nevertheless, substantial amount of computer time is required 
for evaluating the mode sums given in (7.54)-(7.57) making the solution impractical 
unless the convergence of the sums is improved. One way to achieve this is by using 
roof-top basis functions considered next. 


Solution with Piecewise Linear Basis Functions 

The equivalent magnetic current components are now expanded as 

= zjt,M Ir ,T t (x)P,(y ) 

p=l 9=1 

= tf M m P,(x)T,( y) 

p=l 9=1 


(7.64) 
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where 


T.( 0 = 


lull — (s - 1)A£ < £ < sA£ 

A£ 


(s + 1}A£ — £ sA£ < £ < {s + 1)A£ 


(7.65) 


|£ - »M\ > A£ 


for 


= 1,2 5-1 


(7.66) 


and P(0 is the piecewise constant basis function considered before. 

When these are used in (7.39) with the weighting functions the same as the 
expansion functions, we obtain a system similar to (7.44). The external admittance 
elements are now given by 


nr.) = -ji^r r 

^xxJiiTpq ko y (i _i)A v J(,-1)A» J{p- 1)A* 


(7.67) 


(*• + h) 


Go(t‘, r')dx'dxdy dy 


lY .) _ i***tm) r^xM r 


fP^X 

-l)Ax 


(7.68) 


& 1 

dxdy 


G 0 (r ; T')dx'dxdy dy 


(n-.) 


*;w 


« r^'uy) r A ’ r 

k 0 hi- 1)A* J(q-\)£kv J(i-\)Ax J(p 


/(p— l)Ax 


(7.69) 


a 2 


G 0 (r; v')dx dxdy dy 
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—2 jY 0 rti+i)*v 


r{j+ l)Ay /•(»+ l)Ay /■•Ar /-p4 

X.», wL^L 


pAx 

(t — l)Ax V(p— ljAr 



G r 0 (r; r')dx'dxdy'dy 


(7.70) 


The calculation of these elements may be simplified by applying integration by 
parts and sampling the ‘field’ integrands at two points [50]. Taking for example the 
first integral, we have after integrating by parts 



^ir l L {« 1. ™ L 

l T!(x) l Ty)G,ix'ix} dy'dy 


(7.71) 


where T a denotes the derivative of T,. Performing midpoint integration 
for the unprimed integrals by sampling at the two points 
(j ~ 2 ) Ay] an< l repeating the process for the other elements yields 


^xxijpq jY 0 k 0 AxAy j(£ + -)J(| + C) + J(f,0 — (f — ^)/(f ~ 1, C) 


+ 


1 

Ax 


[/*(£- 1 , 0 -/*(£ + 1 , 0 ] 


(7.72) 


+ [/({ + 1, C) - 2/«, C) + /(( - 1, 01 J 

_2 jY 

Y*vijpq = ^ M(00 + Ht + 1»0 + ^(0C - 1 ) - 1(( + 1,C — 1)] (7.73) 

_2 jY 

Y v*n n = ~~k^ [ -/ (00 + A0C + 1) + l{i - 1,0 - J{( - i,C + 1)] (7.74) 

Y wijpg = —jY 0 k 0 AxAy j(C + -)/(0C + 1) + !{(,() - (C - |)/(£,C ~ 1) 
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+ 



(£> C — i) — WU + ^)1 


(7.75) 


+ (1^,5 ^«.C + 1) - 2/«.C) + /(4.C - !)]} 

where £ = i — p, £ = j — ? an d 


m,o = 

/■( C+£)^v 

1 - y====dxdy 

'(4-i)Ax 4 ttvx 2 + y z 

(7.76) 

1,(00 = 

4 

r((+k)*v 
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/ ^ " r~rr — 

r (4-l)Ax 47TVX 3 + y z 

(7.77) 


r(t+\)*v 

/■({+i)A* e _jfco >/ l3+w2 , , 

/ ^7 — - 2 dxdy 

'(S-l)Ax 47T-/X 2 + y 3 

(7.78) 


The above integrals can be evaluated using a four term Taylor series expansion and 


axe given in [50]. 

The evaluation of the internal elements of the admittance matrix is straightfor- 
ward and yields the following expressions 
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(7.80) 
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(7.81) 
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(7.82) 


sin (~qAy^j sin (^7 Ay) 


where p mn and C are given in (7.58) and (7.59). The corresponding excitation ele- 
ments are computed as 


(/“%• = 2H ox e’ ko ™ 9 °l iAx ™+°+U- 1 M*v™+<>)AxAy 
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kg Ax sin 6 0 cos <f> ( 


sine 2 ^ 


k„ Ay sin 0 o sin <f> 0 


(7.83) 


(7.84) 


2 I \ 2 

We observe that the asymptotic behavior of the summation terms in (7.79)-(7.82) is 
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now of the form 


1 

(mn)Vn 2 + tti 2 

and the double sums are therefore expected to converge rapidly. The required number 
of modes for convergence within an acceptable tolerance will, of course, depend on 
the geometry and electrical properties of the cavity. 


7.2.3 Results and Validation 

The implementation task of the presented numerical solution is a tedious one as 
is usually the case with most three-dimensional numerical solutions. The validation 
of the code also proved challenging because of the scarcity of reference data and 
the long execution times. The calculation of mode sums constitutes the major part 
of the computer processing time. As noted earlier, for a piecewise constant basis 
implementation the mode sums are slowly converging and Figures 7.2 a and 7.2 b 
give the convergence of the like- polarized and cross-polarized admittance elements, 
respectively, for these basis functions. The double sums were computed using the 
scheme discussed in [85] and the shown curves correspond to the element located at 
the center of a square 1A x 1A aperture. We observe that the mode sums for the 
cross-polarization admittance elements converge rather rapidly. As expected, Y vy of 
the self-cell has not converged even after adding 1000 modes in each direction in 
the piecewise constant basis solution whereas only 50 modes (in each direction) are 
sufficient to reach convergence when using roof-top basis functions as demonstrated 
in Figures 7.2 c and 7.2 d. In general, though, for narrow and long cavities only 
a few modes need be kept along the narrow dimension and this leads to a much 
more rapid convergence since the double sums are essentially reduced to single sums. 
In fact, for very narrow cavities, one may only keep the lowest order mode [86, 
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87]. This was explored to some extent in Chapter 6. It should be noted, though, 
that for long and narrow apertures the pulse basis formulation is preferable to the 
roof-top one described here unless the external self-cell admittance elements are 
more accurately evaluated (i.e., midpoint integration should be replaced with a more 
accurate integration scheme.) 

To validate the presented moment method full-wave formulations and associated 
computer codes we relied on comparisons with data obtained from a corresponding 
finite element-boundary integral (FEM) solution [88]. This was developed in par- 
allel with the moment method/modal solution in an effort to avoid cavity shape 
restrictions and the long processing time required for filling the MoM matrix. In 
the Figures to follow, the RCS pattern is presented for the principal plane cuts of 
the cavity-backed aperture. Figure 7.3 presents the two like- and cross-polarization 
backscatter RCS patterns for a 1.73A deep cavity with a = OTA and b = 0.1A. These 
are conical cuts and were generated with the code based on the piecewise constant 
basis formulation. They are clearly in good agreement with the FEM data and have 

also been found to agree with the only other [67] available calculations that appeared 
recently in the literature. 

Backscatter curves for a filled cavity are given in Figure 7.4. These were generated 
with the code based on the roof-top basis formulation and correspond to a 0.4A x 
0.4A cavity backed aperture, 0.25A deep and filled with homogeneous material having 
e r = 2 — j0.5 and p r = 1.2-jO.l. The principal plane like-polarized RCS patterns 
are again in good agreement with FEM data. Additional curves for a long and 
narrow 2.5A x 0.25A cavity are given in Figures 7.5 and 7.6. They are based on the 
piecewise constant basis formulation and correspond to a 0.25A deep cavity, empty or 
filled with material having e r = 7 - jl.5 and = 1.8-jO.l. Figure 7.5 presents the 
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like-polarized backscatter RCS patterns for the empty cavity in both principal planes 
whereas Figure 7.6 includes the corresponding patterns for the filled cavity. As seen, 
the RCS patterns in the principal plane normal to the long side agree with the scaled 
two-dimensional RCS data (using the conversion (1.26)), whereas the principal plane 
patterns normal to the short side agree with the FEM data. Finally, the scattering 
characteristics of a square cavity (1A x 1A x 0.5A) filled with a high contrast material 
(e r = 12 — y 2.5, fi r = 4.5 - j 1.2) is shown in Figure 7.7. It is noted that for this 
particular case, a sampling interval of A/15 and a total of only 50 modes in each 
directions were sufficient for the MoM solution to reach the converged solution. This 
is, of course, due to the fact that higher order modes are suppressed because of the 
high losses in the material filling the cavity. 
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Figure 7.3: Comparison of conical (6 = 40°) backscatter RCS patterns for a 0.7 A x 
0.1 A x 1.73A empty cavity obtained from the moment method solution 
using piecewise constant basis functions and the finite element method 
(FEM) [67]. 
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Scattering from a Filled Cavity 

a=0.4A , b=0.4A , c=0J25A^ fc =2-j0.5, 



Angle of Incidence 6, deg. 


Figure 7.4: Backscatter RCS elevation patterns for a 0.4A x 0.4 A x 0.25A cavity filled 
with a homogeneous material (c, = 2-j0.5,p r = 1.2-j0.1); Comparison 
of the MoM solution using piecewise linear (roof-top) basis functions with 
the FEM [67]. 
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Figure 7.5: Backscatter RCS elevation patterns for a 2.5A x 2.5A x 0.25A empty 
cavity using piecewise constant basis functions, (a) <f> = <f> 0 = 0 (symbols 
denote FEM results [88]). (b) <f> = <f> 0 = tt/2 (symbols denote the scaled 
two-dimensional RCS data). 
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Figure 7.6: Backscatter RCS elevation patterns for a 2.5A x 2.5A x 0.25 A filled cavity 
( £r = 7 _ jl.5, (i T = 1.8 — ;0.1) using piecewise constant basis functions, 
(a) 4> = <j> 0 = 0 (symbols denote FEM results [88]). (b) <f> = <j> 0 = */ 2 
(symbols denote the scaled two-dimensional RCS data). 
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Scattering from a Filled Cavity 

a=lA , b-lA , c=0.5A 



Figure 7.7: Backscatter RCS elevation patterns for a 1A x 1A x 0.5A cavity filled with 
a high contrast material (e r = 12 - j2.5,/* r = 4.5 - jl.2); Comparison 
of the MoM solution using piecewise constant basis functions with the 

FEM [88]. 
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7.3 GIBC Formulation 

The full-wave formulation presented in the previous sections is applicable only 
for the cavities of rectangular shapes and may not be employed for the analysis of 
nonrectangular geometries. Furthermore, this formulation is not very efficient when 
considering aperture areas larger than 1A 2 . Therefore, in this section, as in the two- 
dimensional case, we present an approximate formulation based on a simulation of 
the cavity backed apertures by an impedance insert satisfying impedance boundary 
conditions. The analysis of three-dimensional impedance inserts was given in Chapter 
5. 

Considering first the SIBC, the relevant integral equations are (5.59) and (5.60). 
Again, by setting F 2 = 0 in (5.64) and (5.65), we have 

F 1 (x,y)M x (x,y) + ^DFT" 1 % - D 2 X ) - M V D X D V ] (} 

= 2Z o H ox exp{jk 0 [sin 0 O (:rcos fa + ysin 4> 0 )]} 

(7.85) 

and 

F l [x,y)M,(x,y) + ^DFT- { [-M.D.D, 4 - D])] (} 


= 2Z 0 Hoy exp {jk a [sin 0 O ( x cos <f> 0 + y sin 4>o )\ } 

(7.86) 

which may be solved by the CGFFT method. 

Figure 7.8 shows a comparison of the full- wave and SIBC solutions for scattering 
from the filled 1A 2 square cavity considered earlier. Since the filling material is of 
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high contrast (e r = 12-j2.5,^ = 4.5 1.2), good agreement is observed between 

the two solutions. When the losses in the cavity are not sufficiently high, the SIBC 
is not applicable and a higher order GIBC is required. Figure 7.9 shows the results 
for a long cavity 2.5A x 2.5A x 0.25A filled by a material of lower loss (e r = 7 - jl.5 
and n T = 1.8 - jO.l). The second order GIBC based on equations (5.64) and (5.65) 
is seen to yield an improved result for the longitudinal cut considered in this case. 
For H-polarization incidence, the magnetic currents do not vanish at i = ±1.25 A 
and thus, as noted for two-dimensional cavities, the GIBC simulation would not be 
of acceptable accuracy in predicting the currents in the immediate vicinity of the 
cavity edges. This is a limitation in the numerical and analytical application of 
the GIBC, and stems from their non-uniqueness [89]. As noted in [90], additional 
conditions must be imposed at the cavity terminations to supplement the GIBC. 
Although the notion of these supplemental conditions is understood, their numerical 
im plementation is cumbersome and inefficient in the context of the CGFFT solution. 
Thus, additional research is required before the GIBC can be employed for simulating 
coatings and filled cavities with abrupt terminations. On the other hand, if the cavity 
depth on the coating thickness is tapered to zero-as is often the case in practice-the 
presented formulation is then directly applicable. Unfortunately, no reference data 
are available for tapered three-dimensional coatings and cavities which will permit 
validation of the GIBC formulation. 

7.4 Summary 

A full- wave moment method formulation was presented for computing the scatter- 
ing by an aperture formed by a rectangular cavity in a ground plane. In constructing 
the integral equations, the equivalence principle was employed to introduce equiva- 
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Scattering from a Filled Cavity 

a=lA , b=lA , c=0.5A 



Figure 7.8: Backscatter RCS elevation patterns for a 1A x 1A x 0.5A cavity filled with 
a high contrast material (e, = 12 - j2.5,/z r = 4.5 — >1.2); Comparison of 
the full-wave (MoM) with the SIBC (CGFFT) solutions 
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Figure 7 9: Comparison of the E-polarization scattering patterns for the long 2.5 A x 
2.5A x 0.25A filled cavity (e, = 7 - j 1.5, ** = 1.8 - j'0.1) as obtained by 
the full- wave and approximate formulations (longitudinal 4> — 0° cut). 
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lent magnetic currents across the cavity aperture. The fields interior and exterior to 
the cavity were then expressed as the radiation of the equivalent magnetic currents 
in conjunction with the modal and free space Green’s function, respectively. Cou- 
pled integral equations for the two components of the magnetic currents were then 
obtained by enforcing continuity of the tangential fields across the cavity. These 
were discretized using Galerkin’s technique in conjunction with piecewise constant 
and roof-top basis expansion functions. The resulting matrix system was solved by 
LU decomposition. 

The most challenging aspect of the implementation was the computation of the 
slowly converging mode sums required for the evaluation of the interior admittance 
elements. This is particularly so for the piecewise constant basis implementation 
unless the cavity is narrow in one direction. The roof-top basis implementation pro- 
vided a much more rapid convergence at the expense of complexity in the evaluation 
of the external admittance elements. Nevertheless, as is usually the case with three- 
dimensional moment method solutions, the presented solution demands excessive 
CPU time when the aperture size is beyond one square wavelength. It is, therefore, 

more applicable for smaller cavities and particularly those which are narrow in one 
direction. 

Next, an alternative approach based on the GIBC was presented. In this case, the 
structure is essentially modelled as an impedance insert find since this formulation 
is amenable to a CGFFT solution, larger cavity-backed apertures can be handled. 
For the analysis of cavities filled with electrically dense materials, very good results 
were obtained when employing the SIBC. A second order GIBC was also considered 
which is believed to yield acceptable results for lower contrast material fillings. 

An important contribution of this Chapter was the presentation of RCS patterns 
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for various empty and filled rectangular cavities and to our knowledge these are the 
first validated patterns to appear in the literature for this basic cavity shape. 
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Part III 

VECTOR-CONCURRENT 

APPLICATIONS 



CHAPTER VIII 


OPTIMIZATION OF THE CGFFT 
ALGORITHM 


8.1 Introduction 

Computational electromagnetics relies heavily on vector-oriented algorithms to 
simulate complex problems. With the computer technology approaching the limits 
of semiconductor speeds, the exploitation of parallel processing has emerged in order 
to meet the processing demands of computationally intensive applications in electro- 
magnetics. Most modern computing facilities now offer vector and parallel processing 
capabilities. A vector facility exploits the independence of operations, particularly 
those associated with the elements in an array or vector. In such machines, instruc- 
tions are vectorized and distributed across different vector processors for concurrent 
execution as opposed to the traditional approach where the computers are limited 
to sequential processing of data on a single scalar processing unit. 

The CGFFT lends itself to efficient execution in vectorized fashion. Most oper- 
ations involve array manipulations which are vectorizable. Also, several of the steps 
in the iteration algorithm can be treated independently and can thus be performed 
on different processors. Most importantly, since the FFT is a highly vectorizable 
algorithm, it plays a major role in the speed of the solution algorithm and overall 
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efficiency of the optimized code. In this chapter, a vector- concurrent form of the 
CGFFT method suitable for implementation on parallel multiprocessor systems is 
applied to the problem of scattering from electrically large planar structures. To 
demonstrate the speed advantage which can be realized when executing the CGFFT 
solution on a vector-concurrent facility, a few tests were performed on the supercom- 
puters and mini-supercomputers. 

8.2 Optimization 


Before an assessment of the vectorizability of the CGFFT algorithm, a brief 
review of some general concepts in vector and parallel processing will be presented. 
This discussion is followed by an overview of the optimized CGFFT algorithm used 
in this study. 

8.2.1 Vectorization 

The crux of parallel computing is the process of vectorization and distribution 
of code among multiple processors. Typically, a vector instruction is capable of 
operating on 32 to 128 elements of data at once, depending on the machine used, 
resulting in two to four times gain in speed over the corresponding sequential scalar 
instruction. Vectorization requires some degree of independence in the access of 
data by the code. A dependence occurs when two statements— or iterations of the 
same statement— refer to the same storage location. Some data dependences inhibit 
vectorization; they are called recurrences. By changing the structure of the code, it 
may be possible to eliminate a recurrence and vectorize the modified code. 

A typical example of vectorization occurs when performing element by element 
addition or product of two independent arrays/ vectors. In a scalar machine, each 
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element product or addition will be done sequentially, whereas m a vector facility 
vector registers are employed to perform several of the element operations concur- 
rently. That is, when a DO loop is encountered, the loop iterations are not executed 
sequentially but in parallel, provided there are no data dependences among the loop 
iterations. When a parallel (concurrent) facility is also available, independent op- 
erations or sections of the program may be executed on different processors. In 
this manner, several matrix operations involving independent vectors /arrays may be 
performed in parallel. 

In order to measure the improvement in the speed of a vectorized code, several 
parameters are defined. The program speedup is defined as the increase in the speed 
of execution when a code is run in vector mode relative to that in the scalar mode. 
Therefore, referring to Figure 8.1, if a code runs for T. seconds in scalar mode and 
for T v seconds in vector mode (i.e. after optimization), the corresponding program 
speedup is given by 

program speedup = T t /T v (8-1) 

Also, the vector content of a program is that percent of the scalar code which 
vectorizes. Thus, if for a given code, the scalar portion which may not be vectorized 
runs in t, seconds, and that which is vectorizable runs in t v seconds in scalar mode 
and in t v f seconds in vector mode 1 , we have 

vector content = t v jT a x 100 (®-2) 


vector speedup = t v /tvj 

l t v] would be the time the code actually spends in the vector facilities. 
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scalar 

execution, Ts 


vector 

execution, Tv 


Figure 8.1: Scalar and vector execution times in a typical vectorized code. 

Typically, programs with more than about 70% vector content run 1.5 to 2 times 
faster in vector mode. However, as will be shown later, higher vector contents are 
achievable for the case of CGFFT algorithm. 

8.2.2 Concurrency 

Although it is increasingly expensive to make a single processor faster, fairly fast 
processors are inexpensive. Therefore using several relatively inexpensive processors 
in parallel is often more efficient than using a single fast processor. A multiple 
processor system can devote several processors to the execution of different parts 
of a single program simultaneously. The compiler inserts protective synchronization 
code into the optimized loops so that the multiple processors work together without 
interfering with each other. Synchronization is needed to prevent conflicts in the 
use of memory shared by parallel tasks and is considered the major cost of parallel 
processing. Optimization is suppressed whenever the possibility of a data dependence 
exists. Also, from Amdahl s law the increase in the speed of execution in concurrent 
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mode as the number of processors is increased, approaches an upper limit set by the 
presence of the sequential constructs in the algorithm. 

In general, the concepts of data dependence in parallel processing are the same 
as those in vector processing; the consequences of certain dependences, however, are 
different. 

A parameter of interest when processing in concurrent mode is the efficiency 
of execution. Efficiency is a measure of parallelism in the algorithm. An efficient 
parallel tasking system makes possible a nearly linear speedup m performance as 
processors are added, if the algorithm is parallel. Thus, if T (n) denotes the time 
required to run on n processors, we may define 

J ’(1) 

concurrency speedup = (8-4) 

and 

efficiency% = ~j^i) x 1^0 (8-^) 

8.3 Optimized CGFFT Algorithm 

Here, the vectorizable nature of the CGFFT algorithm is exploited by identifying 
the major processes involved in a given iteration. From (2.1), it can be seen that 
each iteration in (2.3) requires two convolution operations, two norm calculations 
and three scalar products. A considerable amount of computation time is spent in 
the calculation of the convolutions carried out in .A[P„] and ,A 0 [R„]. Each convo- 
lution includes a pair of forward and inverse Fourier transform operations on the 
relevant components of the current density vector along with a Hadamard (element 
by element) multiplication of the current and the dyadic Green’s function in the 
spectral domain. Since for a given current component there is no data recurrence 
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at a particular point in an iteration, these operations may be vectorized to increase 
the speed of calculations. The same observation is true for the computation of the 
norms and the dot products used to update the current J, the residual vector R, 
and the search vector P. More importantly, since the FFT is a highly vectorizable 
algorithm, it plays a major role in the speed and efficiency of the optimized code. 

The processing was carried out on two vector machines available at the time of 
this study, namely, the Affiant FX/8 multiprocessor and the IBM 3090/600E super- 
computer. Some general guidelines for code optimization are given in Appendix E 
and have been followed in optimizing the CGFFT algorithm in the present study. 
It should be noted that the data reported here on the actual performance of the 
algorithm on a given mode of execution will be of little value, as current and future 
advancements in computer technology renders them obsolete; however, it is the rel- 
ative performance improvement which is of interest when the algorithm is executed 
in the optimized mode as compared with the sequential mode. 

To assess the efficiency of the optimized code, the rectangular plate problem 
discussed in Chapter 4 was examined in some detail. Tables 8.1 and 8.2 show the 
performance of the optimized CGFFT algorithm executed on the Affiant and IBM 
vector facilities for computing the currents on a 21 x 2A conducting plate. The 
plate was assumed to be illuminated at normal incidence by an E-polarization plane 
wave and 63 x 63 unknowns and 128 x 128 FFT pad (order 1) were employed. The 
percent vectorizable code of 97% indicates a highly optimized algorithm resulting 
in a program speedup of more than four times. It is clear from the tables that the 
vectorized FFT is mainly responsible and contributes to the speedup in the execution 
time. The performance point of the IBM supercomputer for this particular case 
is given in Figure 8.2 indicating the remarkably high efficiency of the optimized 
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algorithm. 

In the case of the Alliant, the overall speedup was more than 600 percent per 
iteration. The speedup in execution time is even more impressive when all four 
processors of the Alliant are utilized. As seen from Table (8.4), a speedup of 3.5 was 
achieved when using four concurrent processors at an efficiency of 88%. This implies 
a speedup of more than 20 times per iteration when combined with the data in Table 
(8.2). Again, the improvements in the performance of the algorithm is attributed 
to the vectorized FFT which is the most significant factor in the solution process. 
This is illustrated in Figure 8.3 where the distribution of the CPU time among the 
computationally intensive routines in the scalar and vector modes are shown. 

In order to further evaluate the efficiency of the method when the size of the 
problem grows, the cases of 5A x 5A and 10A x 10A plates were also considered with 
the corresponding results reported in Tables 8.4 and 8.5. The sampling density in 
these cases were 625 unknowns/ A 2 with a FFT pad of order one. Interestingly, similar 
speedups are observed for larger plates indicating that the aforementioned results 
are independent of the cache memory. Figures 8.4 and 8.5 show the components 
of the surface current densities excited on the conducting plates and calculated by 
the CGFFT method in the vector- concurrent mode. It should be noted that the 
calculation of the surface currents associated with the conducting plate in Figure 8.5 
required 125,000 unknowns. This large number of unknowns presents a challenge 
for direct matrix inversion approaches because of their large storage requirement. In 
contrast, the CGFFT solution could be performed on a relatively small computer. 

Finally, the backscattering behavior of an equilateral triangular conducting plate 
of side length 5A at near grazing (conical cut at 0 o = 0 = 80°) is shown in Figure 
8.6. Due to the symmetry of the problem, the full-range data was replicated from 
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that calculated in the range 0 < <f> < 60°. 

8.4 Summary 

It was shown that the conjugate gradient FFT algorithm is suitable for vector- 
concurrent optimization and may be efficiently implemented on multi-processor com- 
puters. The FFT plays a crucial role in the speedup and the efficiency of such an 
application. As the size of the problem becomes larger, there is a corresponding 
degradation in the performance of the optimized code due to the complexity of the 
memory cache references. This complexity, however, has not proven restrictive in 
the examples considered because the CGFFT method does not suffer from the same 
memory requirements as the direct methods do. Thus, relatively large problems can 
be handled without considerable loss of efficiency and speed. 

Although this study was concerned with automatic parallelization, which is lim- 
ited to optimization of individual loops, parallelism at a larger granularity can be 
specified by the programmer to achieve a superior performance for more complex 
problems. An example of such an application is the problem of scattering by a di- 
electric plate of finite thickness where the normal component of the current density 
is totally independent of the planar components and can be solved for by a dedicated 
processor in parallel with them. 



CGFFT Code 


Execution Mode 


Performance 

Scalar 

(scalar FFT) 

Vector 
(scalar FFT) 

Vector 

(vector FFT) 

ELAPSED CPU, sec 

148 

129 

34 

VECTOR CPU, sec 

- 

9 

30 

VECTORIZABLE CODE, sec 

- 

28 

144 

VECTOR CONTENT 

- 

18.9% 

97.3% 

VECTOR SPEEDUP 

- 

3.1 

4.8 

PROGRAM SPEEDUP 

- 

1.15 

4.35 


Table 8.1: Performance of the scalar and vectorized code on the IBM 3090. 


CGFFT Code 


Execution Mode 


Performance 


Scalar 

(scalar FFT) 


Vector 
(scalar FFT) 


Vector 

(vector FFT) 


INITIALIZATION, sec 
CGFFT LOOP, sec 
TOTAL CPU TIME, sec 
ITERATIONS 
PER ITERATION, sec 


2.39 

4307.0 

4309.3 

111 

38.80 


0.72 

1255.2 

1255.9 

111 

11.31 


0.80 

342.8 

343.6 

59 

5.8 


MEGAFLOPS 
PROGRAM SPEEDUP 
SPEEDUP /ITERATION 


0.0527 


0.1807 

3.43 

3.43 


0.3512 

12.54 

6.69 


Table 8.2: Performance of the scalar and vectorized code on the Alliant FX/8. 
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Optimized CGFFT 

NO. of processors 

Performance 

i 

2 

3 

4 

INITIALIZATION, sec 

1.60 

0.86 

0.61 

0.50 

CGFFT LOOP, sec 

407.4 

211.0 

146.7 

115.6 

TOTAL CPU TIME, sec 

409.0 

211.8 

147.3 

116.1 

ITERATIONS 

59 

59 

59 

59 

PER ITERATION, sec 

6.91 

3.58 

2.49 

1.96 

MEGAFLOPS 

1 



1.0340 

SPEEDUP 


1.93 

2.78 

3.52 

EFFICIENCY 

■1 

96.5% 

92.7% 

88.0% 


Table 8.3: Vector- Concurrent performance for a 2A x 2A plate. 
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% Vectorized 


Design Point 
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<D 

o> 

a* 

m 


cG 
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bO 

O 

tH 

Ph 



4 5 6 7 8 

Vector Speedup 


Figure 8.2: Performance of the optimized CGFFT algorithm on the IBM 3090. 
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Speed Improvement of the Vectorized Code 



MAIN NORM OPERATOR FFT TOTAL 


Routines 


Figure 8.3: Distribution of the CPU time among the computationally intensive rou- 
tines. 



NO. of processors 


6.35 3.26 2.32 1.90 

1096.6 568.9 398.4 319.0 


Optimized CGFFT 


Performance 


INITIALIZATION, sec 
CGFFT LOOP, sec 
TOTAL CPU TIME, sec 
ITERATIONS 
PER ITERATION, sec 


MEGAFLOPS 

SPEEDUP 

EFFICIENCY 


Table 8.4: Vector-Concurrent performance for a 5A x 5A plate 


1103.0 

46 

23.98 


0.3848 


400.7 320.9 

46 46 

8.71 6.98 


1.0590 

1.3226 

2.75 

3.43 

91.7% 

85.8% 


Optimized 

NO. of processors 

Performance 

1 

2 

3 

* 1 

INITIALIZATION, sec 

24.55 

13.60 

9.10 

7.23 

CGFFT LOOP, sec 

3673.3 

1973.0 

1349.8 

1080.5 

TOTAL CPU TIME, sec 

3697.8 

1986.6 

1358.9 

1087.7 

ITERATIONS 

38 

38 

38 

38 

PER ITERATION, sec 

97.31 

52.28 

35.76 

28.62 


MEGAFLOPS 0.4223 

SPEEDUP 

EFFICIENCY - 


Table 8.5: Vector- Concurrent performance for a 10A x 10A plate 















Figure 8.4: E- polarization plane wave scattering from a 5A x 5A conducting plate at 
normal incidence (125 x 125 unknowns and FFT pad of order g = 1). 
(a) Co-polarized component of the current density, (b) Cross-polarized 
component of the current density. 



(b) 
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Figure 8.5: E-polarization plane wave scattering from a 10A x 10A conducting plate 
at normal incidence (250 x 250 unknowns and FFT pad of order Q = 1). 

(a) Co-polarized component of the current density, (b) Cross- polarized 
component of the current density. 



(b) 
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Scattering from a 5A Equilateral Triangular Plate 
Conical Cut; 0=80°. 



Figure 8.6: Conical backscattering cross section of an equilateral triangular conduct- 
ing plate of side 5A at 6 = 80°; 625 unknowns/A 2 ; Maximum number of 
iterations: 250; Average tolerance: 0.015 (E-Pol.) and 0.025 (H-Pol.). 



CHAPTER IX 


Conclusions 


The theoretical and computational aspects related to the application of the con- 
jugate gradient FFT method in computational electromagnetics have been examined. 

The first Part of the thesis was devoted to the problems of electromagnetic radi- 
ation and scattering from linear, cylindrical, and planar structures. Both perfectly 
conducting and imperfect bodies were treated. A number of highly efficient and 
accurate numerical codes have been developed for the solution of these problems. 
These programs cover a broad range of operations in terms of frequency, material 
composition, and structural geometry and may be used for both analysis and design 

purposes. 

The provisions of incorporating various expansion functions into the CGFFT 
method was discussed in Chapter 2. It was found that by employing subdomain 
basis functions, the convergence rate of the CGFFT method can be improved drasti- 
cally. In particular, a quantitative measure of convergence improvement was estab- 
lished for a class of such basis functions. Illustrative examples of CGFFT applica- 
tions to two- and three-dimensional problems were presented in Chapters 3 and 4, 
respectively. Two different but related approaches in computing the mtegrodiffer- 
ential convolutions were presented. The first approach (CGFT) was based on em- 
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ploying the sampled continuous transform of the Green’s function, while the other 
(CGDFT) employed finite duration discrete Fourier transforms. The latter method 
was found to provide a more efficient (and more accurate) simulation, particularly 
for E-polarization cases. 

Other variations in formulating the CGFFT method have appeared in the liter- 
ature recently. These approaches differ primarily in the chosen (or implied) basis 
functions for the unknown current density, and the method of computing the inte- 
grodifferential operator in the spectral domain. Some examples of these variations 
have already been discussed in this thesis. In addition to these approaches, some 
authors have proposed other methods of formulating the integral equations. In par- 
ticular, it has been shown that by introducing the surface charge density in the 
integral formulation and expanding the vector potential instead of the current den- 
sity, smoother and more stable solutions may be achieved [91]. 

Incorporation of the impedance boundary conditions into CGFFT was consid- 
ered in Chapter 5. Advantages include the elimination of a need to sample within 
the volume, further reducing the memory demand and also avoiding the difficulties 
associated with the calculation of the Green’s function. Application of the gener- 
alized impedance boundary conditions in modelling the two- and three-dimensional 
impedance inserts was shown to be compatible with the basic CGFFT formulation, 

thus allowing an efficient simulation of coated structures and filled cavity-backed 
apertures. 

The problems of scattering from two- and three-dimensional rectangular grooves 
and cavities were studied in Part Two of the thesis. After presenting a general full- 
wave analysis, approximate solutions based on the impedance boundary conditions 
were considered. GIBCs of up to order 3 were employed. The formulations based 
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on these boundary conditions were amenable to a CGFFT solution and were found 
easier to implement than the full-wave formulation. The predicted currents based 
on the GIBC simulations are in general incorrect near the edges. However, for high- 
contrast material fillings, the SIBC was found adequate in modeling the groove. 
Further research is needed to study the applicability of the higher order conditions 
to terminated cavity and coated structures, and to establish the ranges of validity for 
such applications. For the two-dimensional case, a hybrid exact-GIBC approach was 
proposed which provided a good prediction of the scattering behavior of rectangular 
grooves without compromising the advantages offered by CGFFT. In contrast, when 
the impedance variations of the insert are sufficiently slow, a direct application of 

the impedance boundary condition is sufficient. 

A vector-concurrent implementation of the CGFFT method was presented in Part 
Three. It was shown that the CGFFT algorithm is highly vectorizable and may be 
efficiently implemented on supercomputers and multiprocessor machines. Vectoriza- 
tion and parallelization of the underlying algorithms will be of great importance in 
reducing the computation time and improving the efficiency of the CGFFT solution 

method. 

Beyond the basic applications considered in this study, the extension to three- 
dimensional structures with and without anisotropy is straightforward by employing 
the three-dimensional FFT. Also, the CGFFT method is directly applicable in solving 
systems relating to scattering, transmission and radiation by periodic structures and 
arrays. In that case, the resulting system is discrete and no need arises for corrective 
measures due to discretization. 

New emerging methodologies which combine the CGFFT with other numerical 
techniques are of potential importance in future research. In addition to the hy- 
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bnd exact- GIBC solution discussed in Chapter 6, the Finite Element and CGFFT 
methods can be combined to reduce the dimensionality of the required FFT and 
consequently improve the efficiency of the solution process [92]. 

The numerical results presented in this thesis will serve two purposes. First, they 
can be used as reference for future developments in this area. Second, future work 
may use the various programs developed in this study to investigate the behavior 
of different scatterers in an attempt to develop simple mathematical and physical 
models. Since the main advantage of the CGFFT method in comparison with matrix 
inversion techniques is its reduced memory demand, it is particularly useful in large 
scale electromagnetic simulations. 
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APPENDIX A 

THE METHOD OF MOMENTS 

Traditionally, equation ( 1.1) is solved directly by the Method of Moments(MoM) [20, 
93]. The method of moments is a projective method in which a functional equation 
in an infinite dimensional function space is approximated by a matrix equation in a 
finite dimensional subspace. 

Consider the linear operator equation 

A [f) = 9 (A.l) 

where A is the linear operator, g is a known function, and / is an unknown function 
to be determined. In the method of moments the unknown function / is represented 

approximately by a linear combination of a finite set of functions f n in the domain 
of A 

N 

f - H C n/n fn € V A (A.2) 

n=l 

where c„ are scalars to be determined. The functions /„ are known as basis or 
expansion functions. Substituting (A.2) into (A.l) gives 

N 

H Cn^[/n] ^ 9 (A. 3) 

n=l 

where the linearity of the operator has been employed. Defining the residual error R 

R = 9~Y. c./t[/„| 

n=l 


(A.4) 
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the coefficients c„ are computed so that the residual error is orthogonalized, with 
respect to an inner product, to a sequence of weighting functions w n defined in the 

range of A 

< w m ,R >= 0 m = w m eH A (A.5) 

The above inner products are called the weighted residuals. This represents a system 

of linear equations 

N 

y: Cti < W m ,A[fn} > — < W m ,g > 

n= 1 

which can also be put in the matrix form 


m 




(A.6) 


[^mn][^n) "" [ 9m\ 


(A.7) 


where [A mn ] is the matrix of elements 


A mn =< l V m ,A[f n ] > 


(A. 8) 


and [^ m ] is the column vector 

g m =< u>m,g > 

If [A mn ] is nonsingular, its inverse exists, and [c„] is given by 

[fin] = [A m n] [5m] 

The solution / is then obtained from (A. 2) 

/ = [/n] < [Amn] 1 [5m] 


(A.9) 


(A.10) 


(A. 11) 


where [f n Y is the row vector of basis functions. 

Two classical approaches have found utility in choosing the weighting functions 
te n . These are referred to as the Galerkin’s method and the point matching method. 
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A.l Galerkin’s Method 

Galerkin’s method may be considered as the specialization of moment method to 

the case of self-adjoint operators. The adjoint operator A a is defined with respect to 
the inner product as 

< w, A[f ] >=< A a [w), f> feV A weV a A (A. 12) 

and if the domains of A and A a are the same, we can choose w n = /„. For self-adjoint 
operators (A = A a ), this choice of weighting functions makes [A mn ] a symmetric 
matrix which may be desirable from a numerical standpoint. 

A. 2 Point Matching 

If the weighting functions are formally chosen to be Dirac delta functions, equa- 
tion (A. 3) is satisfied at discrete points in the region of interest. This is the simplest 
specialization of the moment method. The major advantage of this method is that 
the integrations represented by the inner products (A.8) and (A. 9) now become triv- 
ial since they are evaluated at discrete points. 
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appendix b 

THE FREE SPACE GREEN’S FUNCTION 
AND ITS TRANSFORM 


Consider the complex vector wave equation (Helmholtz equation) satisfied by the 
electric field in a homogeneous isotropic medium 

V x V x E - k 2 E = -jw/iJ ( B>1 ) 

The field may be expressed in terms of the Hertz vector potential of electric type 

n(r) = -ij M v 3«)G(rSW < B ' 2 ' 

where G is the scalar free space Green’s function satisfying the scalar wave equation 

V 2 G(r) + k 2 G{r) = -£(r) ( B - 3 ) 

The Green’s function can be regarded as the response due to a point source and it 
is of interest to find the Green’s function and its Fourier transform corresponding to 
a line source (two dimensional case) and a point source (three dimensional case). 

B.l Line Source 

For a two dimensional problem (|j = 0), the Helmholtz equation reads 

(Jj; + ^5 + k 2 )G(i,y) = -S(x)S(y) 


(B.4) 
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Due to the axial symmetry of the problem, the wave equation can be written 
cylindrical coordinates as 


m 


^dp 2 + pdp + ~ 


(B.5) 


where 6(p) = 6(x)6(y). Outside the source region, the right hand side is zero and 
(B.5) is the Bessel equation of the zeroth order. Therefore, in view of the time 
convention e 3U,t , a solution of (B.5) representing an outgoing wave that satisfies the 
radiation condition is the Hankel function of the second kind [<X>f 


G(p) = 

V 


aS 


Xhe asymptotic expansion of the Hankel function is given by [44] 


~ kp - 

and G shows the proper behavior in the far field. 

In order to find the Fourier transform G, we write 


oo 


G(x ’ y) ~ (2^ C L G{k " 

6(x ’ v) - ?2 

and substitute these in (B.4) to obtain for all x and y we have 

7SiCC eiit ‘’* k "‘' dk ^’ Vx ’* 


(2tt) 


(2ic)- j — oo J — oo 

Consequently, we may formally write 


(B.6) 


(B.7) 


(B.8) 

(B.9) 


(B.10) 


G(k x , k y ) = 


V-k 2 -k 2 


(B.ll) 
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and equivalently, 


r , t A l_ r r — - — ^ x ^ y) dk s dk v (B.i2) 

G(x,y) - ~ (2ir )2 L" 1,0 V - k* - & 


(2tt) 

However, this integral is undefined for real values of k, because the poles 

k y = ±y[k* - 


(B.13) 


are located in the path of integration on the real fc„-axis (Figure B.l). This difficulty 
can be alleviated by introducing a small loss in the medium so that k = k' - jk", 
and the poles are given by 

k v = T(k' p -jty (B.14) 


The singularities are now located off the real axis and the inverse transform (B.12) 
is defined. Furthermore, in accordance with the radiation condition, we demand the 
following functional form for the transform 

G+(k„y) = K>0 (B.15) 

G-(k„y) = r(k,)e’ k "\ y<0 ( B ' 16 > 


where fc* are such that 

*e{fc y + } < 0 

i 

5m{A:jJ' } > 0 


and 


*e{*:;} > 0 
3m{A:~} < 0 


(B.17) 


The k y integral may now be carried out by contour integration in the upper and 
lower half planes corresponding to y > 0 and y < 0, respectively. Using Jordan s 
Lemma and Cauchy’s theorem, we obtain 


G(x,y) = 

1 

f oo e i*»v .. 

Lhf‘ dk - 

o 

A 

4 7TJ j 

G{x,y) = 

1 

4x7 - 

r -~~e^ klX dk x , 

1-00 Ky 

A 

o 


(B.19) 
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Figure B.l: Path of integration for the Fourier transform integral, 
where axe given by (B.13): 

k t = (B.20) 

Therefore, by uniqueness of the solution to partial differential equations satisfying 
the required boundary conditions, the final result holding for all values of y is given 

by 


j-H^(kp) = -i- r° - - e jkxX dk 

K 4nj J-oo sen (v)k.. C 1 


oo g-ikirlvl 


4 irj J — OO 


Vy 


(B.21) 


where sgn is the signum function and k v satisfies the second of conditions (B.17). 

The above equation was derived for a lossy medium ( k complex) but it remains 
valid for the lossless case (Ar real) provided the path of integration is stipulated to 
avoid the singularities [94]. In this case k v should satisfy 



yfi* - k l k > k x 
~jy/ k l - k 2 k <k x 


(B.22) 



250 


In the limit as \y\ -* 0, we have 


i ■p 1 | 3? e {^y} > ® 

4 J 2jky I < 0 


(B.23) 


and when k is real, the Fourier transform pair is given by 

1 

I —r 


G(fc|x|) 


B.2 Point Source 


k > k x 


2jyjk*-kl 

| — 1 - k < k x 

l 2sjkl - V 


(B.24) 


In the three dimensionai case, the Helmholtz equation takes the form 

(~ + + Jrs + k2 ) G ( x ^' z ) - -*( x )%K( 2 ) 

y dx 2 dy 2 oy 2 

Solving this equation in the spherical coordinates, gives 

e -jkr 


(B.25) 


G(r) = 


4xr 


(B.26) 


Following an analysis similar to the two dimensional case, the inverse Fourier trans- 
form integral is introduced as 


G(x,y,z) = -phr r r 

v (27r) J J-oo J-oo J-oo 


(B.27) 


«(».».*) = (B - 28) 
which upon substitution in (B.25) yields the formal expression for G 

nt s 1 r r r - e’ [k * x+k * v+k '* ) dk s dk y dk x 

G(x, y, z) = J_ x L" k 2_ k 2_ k 2_ k f 


(B.29) 
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The presence of the poles 


lc z = - kl - kl 


(B.30) 


in the path of integration on the real fc*-axis renders the above integral undefined. 
Again, introducing a small loss in the medium and following the procedure outlined 
in the previous section, the k x integration may be carried out in the complex plane. 
Thus, contour integration in the upper and lower half planes corresponding to z > 0 
and z < 0, respectively yields 


G{ x i Vi z ) 


j r°° y°° 

8tT 2 J-oo J-oo 






■e* k * x + k '*Uk x dk y , 


z > 0 


(B.31) 


G[x ' y ’ z) “ 8 7 Li 

where k*” 1 are given by 


°° ^e j{k * x+k * v) dk x dk v , z < 0 (B.32) 


* t = - kl - k^ y 


(B.33) 


and they satisfy conditions similar to (B.17). Combining the two equations, 
obtain for all z 


we 


e jkr j 2 e -i*,|z| 

4^r “ 2j (2tt) 2 J-oo J-oo sgn {z)k z ej(k ’ r+hvV)dkxdky ' V * ( B - 34 ) 


where 


(B.35) 


3Re{/:*} > 0 
3m{£,} < 0 

When k z is real (lossless case), the path of integration is deformed to exclude the 
real poles and (B.34) remains valid provided 

_/ k’>kl + kl 


k,= 


-jfti + - * ! k*<kl + P 


(B.36) 
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In the limit as \z\ -*■ 0, we have 

e -jks/x 2 + y 2 i 

4tt \/x 2 + ip 2 j k z 


and when k is real the Fourier transform pair is given by 


G{x,y) 


2 jy[& -k 2 - ki 

1 

l 2y/k> + k*-V 


k 2 >kl + k l 

P<kl+ k 2 y 


(B.37) 


(B.38) 
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APPENDIX C 

MATHEMATICAL PROOFS 


C.l Proof of the Transformation (5.10) 

Consider the generalized impedance boundary condition (5.1) applied to the nor- 
mal field component U at the top surface (y = t) of the layer under study 

a m d m U 

f£o (-7*o) m dy m v=t~ 

where a m and a' m are the GIBC constants specific to the material properties of the 

layer. It is desired to replace this condition by an analogous one, applicable at the 
reference plane ( y = 0) as 

v-' -A m d™ U | 

ho ~dir\v=o = ( c - 2 ) 

This may be the case when one is interested in applying the image theory for a 
coating of thickness t on a ground plane. 

In order to find the relation between the two sets of coefficients A m and a m , we 
expand U at y = t and its normal derivatives there in terms of the corresponding 
quantities at y = 0. Hence, by invoking Maclaurin series expansion of f/, we have 



= S(S|^) t,(0, + O ( < " +, ' a *' +, ) (C-3) 
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where we have retained only derivatives of at most M-th order as the original condi- 
tion does not include higher order derivatives. Similarly, we have for the first normal 

derivative 


dUjt) 

dy 


= f 1/(0) + 


and generalizing this to higher orders 

<mt) _ f ( mo) + o(t M+1 - p ,^ +1 ) 

dy p ^\(n-p)!%V 

Substituting the above expansion back into (C.l), we obtain 

M M ( /n-m fjn \ 

sr °m _ y' f _J ) i/(o) = o 

h hu u n - m ) ! d y n ) 

which after a simple re-arrangement of terms reduces to 

1 A (-jM) n STU(0) 


= 0 


S,H*o) m »«o n! m-n 

The above result is in the desired form (C.2) and by inspection, we find 


(C.4) 


(C.5) 


(C.6) 


(C.7) 


A (~jk 0 ty 

“ h n! 

n=0 




m = 0, . . . , M 


C.2 Proof of the Identity (5.52) 


(C.8) 


We would like to show that 

.. dG(x,y,z',x',y 9 ,z') | 

km 3- . n 

x— 0 oz ,2=0 

is a Dirac Delta function. Consider the identity (Appendix B) 

_ LJ— r r (c.9) 

4ir R 2 (2x) 2 J-oo J-<x> k z 


R = ijxi+yi+z 2 


(C.10) 
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implying the Fourier transform relation 


where 


j gj k g 2 

4t R < > 2 k t 


J\j k l + k l - kl , k a < Jkl + kl 


'V-kl-ki , k 0>y /vTv y 


in accordance with the radiation condition. Differentiating with respect to z, 


§-S(K,k v ) = 


whose limit as the sheet is approached is 


lim --e* ktZ = — - 
*- o 2 2 


Upon inverse transformation of the last equation, the desired result is obtained 


lim - ' g ( g » h Z 1 1 u x 

i 1 ® * s =“2 tf(z,y) 


(C.15) 
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appendix d 

OPTIMIZATION TECHNIQUES 


Once a well written code is compiled, most vector and parallel constructs are rec- 
ognized by the compiler and automatically optimized for efficient execution. This is 
achieved through loop level concurrency and vectorization of the sequential code. A 
typical compiler’s optimization strategy may be summarized as follows: 

For each, innermost loop : 

• If vectorizable then 

o If next outer loop is parallelizable then 

★ Concurrent- Outer Vector-Inner 
o Else 

★ Vector- Concurrent 
o Endif 

• Elseif parallelizable then 

★ Scalar-Concurrent 

• Else 

★ Not Optimized 


• Endif 
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The fastest vector- concurrent mode of execution is achieved when all the available 
processors are utilized to attack a single task concurrently with vector operations 
performed on strides of data on each processor. The resulting high-performance, low- 
level parallelism can significantly boost the performance of computationally intensive 
operations such as Fourier transformations. 

In this study the scalar and optimized FFT routines available on the IBM 3090’s 
ESSL library and the Alliant’s Fortran math library were used in both scalar and vec- 
tor modes. These FFT routines are written in conjunction with assembler language 
and generate instructions appropriate for the architecture of the processors-they 
manage data to make efficient uses of the memory hierarchy. 

In addition to automatic optimization, however, most vector compilers provide 
directives for additional control. Compiler directives are user supplied control struc- 
tures to override decisions made by the compiler and to give additional information 
to it. Upon compilation, the directives are interpreted by the processor and con- 
verted to library calls to be executed in a more efficient manner. In particular, the 
associative transformation directive recognizes operations like dot products and norm 

computations as reduction functions and optimizes these otherwise non-vectorizable 
loops. 

Some general guidelines for code optimization in various architectural levels are 
given in Table D.l. 



ARCHITECTURAL PROGRAMMER ACTIONS t 

TECHNIQUES - 


SCALAR PIPELINING Use co mpiler switches (global optimization). 

VECTOR PROCESSING 1. Restructure loops and use compiler directives. 

2. Maximize vector lengths by renesting, merging, 
unrolling loops; largest iterates to be inside. 

3. Eliminate conditionals in loops & distribute them. 

4. “Supervector” loops to fill vector registers. 

5. Move I/O statements out of the loops. 

6. Turn off vectorization for short loops, or some 
vectorized loops with conditio nals, dependences. | 

CONCURRENCY 1. Eliminate/ relocate dependences, scalar 

carry- arounds. 

2. Restructure for Concurrent- Outer Vector- Inner. 

3. Create concurrently-callable subroutines. 

4. Renest /merge loops for more concur rent iterations. 



CACHE/MEMORY 1. Possible problems with strides. 

ACCESS 2. Leftmost array dimension should be the largest. 

4. Outer loop corresponds to the leftmost array index. 

3. Process compact vectors, columns instead of rows. 

4. Localize memory references. 

t Extracted from IBM 3090 and Alliant FX/8 programming manuals 


CACHE/MEMORY 

ACCESS 


Table D.l: Optimization techniques. 
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